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We solve numerically for the first time the two-fluid, Hall-Vinen-Bekarevich-Khalatnikov 
(HVBK) equations for a He-II-like superfluid contained in a differentially rotating, spher- 
ical shell, generalizing previous simulations of viscous spherical Couette flow (SCF) and 
superfluid Taylor-Couette flow. The simulations are conducted for Reynolds numbers in 
the range 1 x 10^ < iZe < 3 x 10**, rotational shear 0.1 < Afi/fi < 0.3, and dimen- 
sionless gap widths 0.2 ^ (5 ^ 0.5. The system tends towards a stationary but unsteady 
state, where the torque oscillates persistently, with amplitude and period determined 
by 6 and ACl/Cl. In axisymnietric superfluid SCF, the number of meridional circulation 
cells multiplies as Re increases, and their shapes become more complex, especially in the 
superfluid component, with multiple secondary cells arising for Re > 10^. The torque 
exerted by the normal component is approximately three times greater in a superfluid 
with anisotropic Hall-Vinen (HV) mutual friction than in a classical viscous fluid or 
a superfluid with isotropic Gortcr-Mellink (CM) mutual friction. HV mutual friction 
also tends to "pinch" meridional circulation cells more than CM mutual friction. The 
boundary condition on the superfluid component, whether no slip or perfect slip, does 
not affect the large-sc;alc structure of the flow appreciably, but it does alter the cores of 
the circulation cells, especially at lower Re. As Re increases, and after initial transients 
die away, the mutual friction force dominates the vortex tension, and the streamlines of 
the superfluid and normal fluid components increasingly resemble each other. In nonax- 
isymmetric superfluid SCF, three-dimensional vortex structures are classifled according 
to topological invariants. For misaligned spheres, the flow is focal throughout most of its 
volume, except for thread-like zones where it is strain-dominated near the equator (invis- 
cid component) and poles (viscous component). A wedge-shaped isosurface of vorticity 
rotates around the equator at roughly the rotation period. For a freely precessing outer 
sphere, the flow is equally strain- and vorticity-dominated throughout its volume. Unsta- 
ble focus/contracting points are slightly more common than stable node/saddle/saddle 
points in the viscous component but not in the inviscid component. Isosurfaces of positive 
and negative vorticity form interlocking poloidal ribbons (viscous component) or toroidal 
tongues (inviscid component) which attach and detach at roughly the rotation period. 
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1. Introduction 



A diverse family of flow states, collectively termed spherical Couette flow (SCF), is 
observed when a viscous fluid fills a differentially rotating, spherical shell. The flow state 
at any instant is determined by the Reynolds number Re, dimensionless gap width S, 
relative angular velocity Afl, and, importantly, the history of the flow. Some of the 
states are steady; others (usually, but not always, those with higher Re, 6, or Ail) are 
unsteady. At low Reynolds numbers {Re < 10'^), the basic flow (0- vortex state) is steady 
and symmetric about the equator. Above a critical Reynolds number, that for small gaps 
{6 < 0.1) can be approximated by Rec ~ 41(1 + 6)6^^^^, a Taylor vortex develops on 
each side of the equator ( Khlebuytin|1968" Junk fc Egbers|2000 1. The meridional velocity 
increases with Re and S { Biihler 1990 Egbers fc Rath 19951, scaling as ve oc Re Ail 
for S < 0.1 and Re 



19861 



For wide gaps (S > 0.3), the flow 
counterrotation; 



< 10^ Yavorskaya et at 
does not develop Taylor vortices except under special conditions [e.g 
see Liu et al. (1996); Loukopoulos fc Karahalios (2004l]. It is unstable with respect to 
nonaxisymmetric perturbations (Belyaev et al. 1978 Yavorskaya et al. 19861. At high 
Reynolds numbers {Re > 10^), the flow develops spiral vortices, shear waves, and her- 
ringbone waves, before entering a fully developed turbulent state as Re increases further 
( [Nakabayashi fc Tsuchida||1988[ |Nakabayashi et aL||2002'^ l. 

The problem of superfluid SCF, for example in He II, has not yet been explored nu- 
merically ( Henderson fc Barenghi|[2004 1 or experimentally. It is not known how the flow 
states differ from viscous SCF, and what transitions are allowed between them. Even in 
cylindrical (Taylor-Couette) geometry, only a limited amount of information exists re- 
garding state transitions in the superfluid problem, for the special cases of very small gaps 
{S ^ 0.02) and small Reynolds numbers {Re < 380) (Henderson et al 



1995 



Henderson 



fc Barenghi 2000 1. Taylor vortices are detected in He 11 at the critical Reynolds num 



bers predicted by linear stability theory {Rec ^ 278) ( Barenghi fc Jones||1988 Barenghi 
19921, but the theoretical predictions are valid only at temperatures T > 2.0 K, close 
to the transition temperature Tc, where the normal fluid component dominates (> 90 
% of the total density). The circulation cells are elongated in the axial direction, and 
anomalous modes (cells rotating in the opposite sense to those in a classical fluid) are 
observed (Henderson fc Barenghi 20001. The streamlines of the normal and superfluid 
components are appreciably different for Re < 10^ but increasingly resemble each other 



as Re increases ( [Henderson fc Barenghi]|1995[ [Peralta et a?.||2006&[ ). 

In this paper, we employ a numerical solver recently developed to solve the Hall-Vinen- 



Bekarevich-Khalatnikov (HVBK) equations for a rotating superfluid ( Peralta et al. 2005 1 



to study the unsteady heh&viom: of SCF in classical (Navier-Stokes) fluids and superfluids, 
in two and three dimensions. First, we perform a set of axisymmetric experiments with 
rotational shear in the range 0.1 ^ Ail/il ^ 0.3 in medium and large gaps (0.2 ^ (5 ^ 0.5). 
The flow is unsteady. The torque, which oscillates persistently and quasiperiodically 
(near but not at the rotation period), can be up to three times greater in a superfluid 
than in a Navier-Stokes fluid at the same Reynolds number. We assemble a partial 



gallery of vortex states, in the same spirit as for classical SCF (Marcus fc Tuckerman 
1987a|& Dumas 1991 Junk fc Egbers 20001; a complete classification lies beyond the 
scope of this paper. Second, we take advantage of the three-dimensional capabilities of 
our numerical solver to investigate two systems that exhibit nonaxisymmetric flow: (i) a 
spherical, differentially rotating shell in which the rotation axes of the inner and outer 
spheres are mutually inclined; and (ii) a spherical, differentially rotating shell in which 
the outer sphere precesses freely, while the inner sphere rotates uniformly or is at rest. 
These systems have never been studied before. We use standard vortex identification 
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methods, introduced by Chong et al. (|1990) in viscous flows, to fully characterize the 



three-dimensional vortex structures we encounter. 

The paper is divided into seven parts. In Section [2] we review the HVBK theory of He 
II. We describe the numerical method in Section [3] and validate it in Section |4l In Section 
|5] we present results for axisymmetric superfluid SCF, empahisizing its time-dependent 
behaviour. We investigate the effects of grid resolution, spectral filtering, superfluid frac- 
tion, (ani)stropic mutual friction, and no-slip/perfect-slip boundary conditions. In Section 
|6] we present results for nonaxisymmetric superfluid SCF for misaligned and precessing 
spheres, emphasizing again the time-dependent behaviour and vortical topology. Labo- 
ratory and astrophysical applications are discussed briefly in Section [7j 



2. HVBK theory 

Hall fc Vinen| ( 1956a|& I and Bekarevich & Khalatnikov ( 1961 1 first devised a two-fluid 



hydrodynamic model to describe rotating He II in the presence of a high density of vortex 
lines with quantized circulation. The HVBK model was rederived by [Hills fc Roberts 



(19771 from first principles, within the framework of classical continuum mechanics. It 



employs thermodynamic variables associated with the fluid as a whole, which satisfy 
conservation equations of mass, momentum, and energy, as in the work of [Green <^ 
Naghdi (19671 and Hills, ( ,1972[ ) on the theory of mixtures. 

In the full HVBK theory, the inertia of the vortex lines is explicitly considered, with 
the superfluid density regarded as an independent thermodynamic variable, resulting in 
a three-fluid set of equations. Vortex line inertia is relevant when studying superfluid 
flow near solid boundaries, as it explicitly includes healing [where the superfluid density 



decreases near a boundary; Donnelly (2005)] and relaxation (which prevents the super- 
fluid fraction from changing instantaneously when the thermodynamic state is altered). 
We do not consider these issues in this paper. Instead, we use the equations of [Hills ^ 
Roberts ( 1977 1 in the HVBK limit where the vortex inertia is zero. 



2.1. HVBK equations of motion 

The incompressible HVBK equations which describe the motion of the superfluid (density 
ps, velocity Vg) and normal fluid (density p„, velocity v„) components take the form ( |Hills| 
fc Roberts[[1977| [Barenghi fc Jones[|l988[ ) 

9v„ 



dt 



+ (v„ • V)v„ = -Vcr„ + J/„V v„ 



-F, 



+ (v, • V)v, = -Vcr, - 



V • v„ = V • V, = 0, 



where (Tc and f7„ are defined as 



as ^ U ~TS ^ ^^(Vn-Vg) H 

p 2p p 



(2.1) 

(2.2) 
(2.3) 

(2.4) 



= c/+^r^+- + ^(v„-v,)^ + 



Here, p is the pressure, p 
Vs is the stiffness parameter (defined in Section 2.2 1, U3 



Ps 



■P 

p^2p' 

is the total density. 



Psl^s\i^s 



(2.5) 



I is the kinematic viscosity, 
V X Vs is the macroscopic 
vorticity (averaged over many vortex lines), and U and S are the internal energy and 
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at a given temperature T. We 



entropy per unit mass, which we take to be uniform 
define the mutual friction F and vortex tension T in (2.1 1 and (2.2| in the next section. 
The incompr essible limit corresponds formally to infinite first and second sound speeds 
Effective pressures Ps and p„ are defined by Vp^ = Vp — 5PnV(v,^ 



( Sonin 
and 



in| |l987| Jt1] 



are defined by Vps — Vp — 5PnV(v^g) 
Vs. In the incompressible limit, only 



the first viscosity coefficient and mutual friction can be included as dissipative processes; 
other transport coefficients involve compression of the normal and superfluid components 
( |Sonin||1987( [Andersson fc Comer||2006| . 



2.2. Mutual friction and vortex tension 

Quantized vortex lines mediate an interaction between the normal fluid and the super- 
fluid component known as mutual friction. The major source of mutual friction in liq- 
uid helium is roton-vortex scattering in the experimentally relevant temperature range 
IK ^ T ^ 2.17 K. For a rectilinear vortex array the mutual friction is anisotropic. 
Hall & Vinen ( 1956a|& I showed experimentally that second sound propagates at differ- 
ent speeds parallel and perpendicular to the rotation axis and is damped in the latter 
direction. They postulated the following form for the mutual friction force per unit mass 
due to a rectilinear vortex array: 



F = ^BlJs X {lJs 



X v„, - 



X 



T). 



(2.6) 



In (2.6), B and B' are temperature-dependent, dimensionless coefficients (Barenghi et al 



19831. The first and second terms on the right-hand side give the force per unit mass 



along and perpendicular to the second sound wave vector respectively. The B coefficient 
attenuates the second sound, while B' shifts its frequency. The term T was not included 
in the original derivation of Hall & Vinen ( 1956a|& I. It was proposed by Andronikashvili 
& Mamaladze (1966), to take into account the curvature of the vortex lines. 

The vortex tension T arises from the local circulation around a vortex line. It was 



added to the HVBK equations of motion by Hills & Roberts (1977) [cf. (2.6)]. Consider a 



vortex line which is slightly curved. The force per unit length, f , which tends to straighten 
the vortex, points towards its centre of curvature and has magnitude e/r, where e is the 
energy per unit length and r is the radius of curvature. In vector form, this can be written 
f — 6(0)5 •V)a>s, with u)s = uJs/\^s\- When extending this argument to many vortex lines, 
the local superfluid velocity around each vortex line is determined by the quantization rule 
f Vs ■ dl — K — h/m, where the integral is calculated around a path enclosing the vortex 
core, m is the mass of the bosonic entity forming the condensate (the helium atom in He 
n, or two neutrons in a neutron superfluid), and h is Planck's constant. The mean area 
density of vortex lines is lOs/k. Hence the average straightening force per unit volume of 
superfluid is {eojs/ k){Cjs ■ V)ws = paVg^^s x (V x lJs), with Vs = e/{psK) ( [Andronikashvili 



& Mamaladze 1966 Khalatnikov 1965). In order to evaluate this force, one needs the 
energy per unit length of vortex line, which is given classically by e = psK^ ln(6o/ao)/(47r), 
where &o is the intervortex spacing and Oq is the core radius of the vortex. The stiffness 
parameter, Vs, in (2.2) is then given by Vs = (K/47r) ln(6o/ao), and the vortex tension 
force per unit mass, T, is written as ( Barenghi fc Jones||1988[ Henderson et a?.||"l995 l 

(2.7) 



Us^^s X (V X Ws). 



f This is a good approximation in neutron stars, for example, an important application where 
the flow is subsonic. Note that we model systems with p„ 7^ 0, which often sustain heat currents. 
However, as long as the flow is slower than the speed of second sound, and no external heat 
source is present, the fluid can be treated as isothermal. 
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Note that has the dimensions of a kinematic viscosity, but its physical meaning is very 
different: it controls the oscillation frequency of Kelvin waves excited on vortex lines 
( [Henderson et a?.||1995[ ). 

Quantized vortices are not always organized into a rectilinear array. If the counterfiow 
speed Vns exceeds a threshold, growing Kelvin waves are excited along the vortex lines 
and the rectilinear array is disrupted to form a self-sustaining, reconnecting, "turbu- 
lent" vortex tangle ( Donnelly|2005[ ) . Experimentally, this is observed in narrow channels 
carrying a heat flux, where second sound waves are attenuated preferentially along v„s 
independently of frequency (and hence velocity gradients) , and the temperature gradient 
is proportional to the cube of the heat flux (Gorter & Mellink 1949 Vinen 1957a|& I. 
These data can be explained by an isotropic mutual friction, called the Gorter-Mellink 
(GM) force. Usually, the GM force per unit volume is written as f = Ap„psU^gV„s, 
where A is a phenomenological constant which is a function of temperature and has val- 
ues 23.1g~^cms ^ A ^ 3310 g~^ cms in the temperature range 1.20 K ^ T ^ 2.16 K 



in liquid helium. Re- writing it as a force per unit mass, as in (2.1 ) and (2.2 1, we have 



F = A' 



PnPsVns 



(2.8) 



where A' = B'^p^Tr^Xi/Sp^xi is a dimensionless, temperature-dependent coefficient, re- 
lated to the original GM constant by A' = Apu, and xi and X2 are dimensionless con- 
stants of order unity ( |Vinen|[r957c Peralta et a?. [2005 1. 



3. Pseudospectral solver 

In this section, we describe our numerical method. We start from a three-dimensional, 
pseudospectral, Navier-Stokes solver, originally developed by S. Balachandar to study 
viscous flows around circular and elliptical cylinders ( Mittal|1995 Mittal & Balachandar 



1995 1, prolate spheroids (Mittal 1999), and rotating spheres (Bagchi & Balachandar 2002 
Giacobello|2005 1. The solver is modified in two steps to solve the Navier-Stokes equation 



in a spherical Couette geometry with time-dependent boundary conditions: 

(a) The absorption filter applied at the outer boundary to enforce outfiow is switched 
off and replaced by a Dirichlet boundary condition (see Section 2.3). The filter smoothly 
attenuates the radial diffusive terms in the Navier-Stokes equations, but it is inappro- 
priate in SCF, which takes place in an enclosed geometry. 

(&) A third-order Adams-Bashforth scheme is used to evolve the fields in time (see 
Section 3.2 1, upgrading the second-order Adams-Bashforth scheme in the original solver. 

The solver is then extended to handle the superfiuid HVBK equations, which are 
mathematically similar to a Navier-Stokes equation coupled to an Euler equation with a 
forcing term. This extension is quite challenging, so we explain the method in enough de- 
tail (in this Section and the Appendices) for the reader to reproduce and verify the results 
if desired. An early attempt to solve the spherical Couette problem with a pseudospectral 
code based on spherical harmonics ( Hollerbach|2000 1 was stymied by numerical instabil- 
ities arising from the sensitivity to boundary conditions Henderson & Barenghi ( 2004 1 , 
R. HoUerbach 2004, private communication]; the basis ftinctions are defined globally, 
so instabilities at the boundaries rapidly contaminate the whole computational domain. 
Our approach, based on restricted Fourier expansions in the angular coordinates and 
Chebyshev polynomials in the radial coordinate, solves these difficulties by combinating 
a low-pass spectral filter ( Don|1994 1 with special boundary conditions for the superfluid 



( [Khalatnikov||1965[ [Hills fc Roberts||1977[ [Peralta et al][2005 1 
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3.1. Geometry 

We consider the motion of an isothermal, incompressible, rotating superfluid, described 



by equations (2.11 and (2.2), contained between two concentric spheres. Points in the 
domain are defined by the spherical coordinates (r, 9, 0), with 



i?2 < ^ i?l, 61 TT, 27r, 



(3.1) 



where Ri and i?2 are the radii of the inner and outer spheres, respectively. The spheres 
are assumed to rotate rigidly, with angular velocities ili{t) and il2{t) respectively. The 
inner sphere rotates about an axis parallel to the z axis; the outer sphere rotates about 
an axis that can be inclined with respect to the z axis, by an angle Oq- The spheres can 
accelerate or decelerate, for example, in response to the back-reaction torque exerted by 



the fluid, or because the outer sphere precesses freely (see Section 3.5.41. All variables 
are made nondimensional using R2 as a unit of length, and as a unit of time, unless 
indicated otherwise. The viscous Reynolds number is defined as i?e = riii?2/j^n and a 
"superfluid" Reynolds number is defined as Rcs = iliR^/i^s- For cases where only the 
inner (outer) sphere rotates, we define Rei — fliRl/vn and Re2 = fl2R2/^n- 

3.2. Algorithm 

The radial coordinate (r) is discretized using a Gauss-Lobatto collocation scheme fBoydJ 



2001 Canuto et al. 1988| ). The angular directions 9 and are discretized uniformly. 
The number of collocation points in the three coordinates is (AV, Ng,N^); their detailed 
coordinate locations are defined in appendix The collocation points are shifted from 
the poles in order to avoid the coordinate singularities at 6* = 0, tt. Note that this 
displacement is small; for a typical grid with Ng — 200, the first grid point is located at 
01 « 7.854 X 10-3 rad. 

In spherical coordinates, the Courant-Friedrichs-Lewy (CFL) stability condition for a 
convective-dominated equation with time step At can be written as a limit on CFL = 
me:>c{ur/ Ar + ug/rA9 + Ucji/r sin 9 Acf)) At, where {Ar, A9, Aip) are the {r,9,(j)) grid spac- 
ings; the maximum is taken over the whole computational domain. Trial and error sug- 
gests that the integrator is stable for CFL < 0.6, as for the Navier-Stokes solver devel- 
oped by Mittal (19991. We usually take 10"^ < Ai ^ 10^"^ in dimensionless units. 

The velocity fields are expanded in terms of Chebyshev polynomials in r and Fourier 
polynomials in 9 and (/). The expansions must obey the pole parity conditions and be 
infinitely differentiable to avoid slow convergence of the numerical scheme and any emer- 



gence of Gibbs phenomena (see Section 3.3 ). The final forms of the expansions (which are 



different for scalars and vectors) are presented in Appendix |B] Differentiation in r and 9 
is performed in physical space, multiplying by a differentiation matrix. Azimuthal deriva- 
tives are calculated in wavenumber space. The explicit form of the r and 9 differentiation 
matrices is presented in Appendix |C.l 



The equations of motion (2.1 1 and (|2.2[) are discretized in time using a time-split algo- 



rithm ( Chorin|1968 Canuto et aL|1988) , using an explicit, third-order, Adams-Bashforth 
method for the non-linear terms and an implicit Crank-Nicolson method for the diffusive 



terms. The final form of the discretized equations is presented in Appendix |C.2| includ- 
ing the pressure correction step. The time-split algorithm is presented in Appendix [C| 
with an explanation of how the original second-order solver is upgraded to third-order 
accuracy. 

3.3. Pole parity 

A number of numerical issues arise when a discretised field is expanded in a Fourier series 



on a sphere ( |Merilees|[T973| |Orszag|[T974l |Mittal||1995"| [Bagchi fc Balachandar]|2002| ) . In 
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particular, the 9 expansion is restricted to the range ^ ^ tt, so a Fourier series 



(periodic in ^ ^ 27r) can only be used with some symmetry restrictions (Bagchi fc 
[Balachandar 20021. In a spherical grid, lines of latitude and longitude intersect at two 



points, and the spherical components of a vector field are discontinuous at the pole even 



when its Cartesian components are continuous ( Swarztrauber||1979 j , so the 9 expansion 
must obey certain boundary conditions at the poles in order to be compatible with the 



expansion. This is called the pole parity problem (Merilees 1973 



Merilees 


1973 


Orszag 


1974 


Yee 



[l981, ,Mittal^l999) . Parity conditions are chosen to ensure that the series expansions are 
difi^erentiable at the poles, avoiding convergence problems and the Gibbs phenomenon 
(Orszag'1974). For a scalar field s{9^ 4>), only certain modes are permitted in the expansion 
( Fornberg 1998 ) . For odd (even) azimuthal wave numbers, the expansion oi s{9, cj)) must 
have odd (even) parity. For a vector field v(0, (p), the r component is continuous across the 
poles, but the 9 and (j) components change sign. The radial component of the vector field 
follows the same parity rule as for a scalar. For the tangential and azimuthal components, 
expansions with odd (even) azimuthal wave numbers must have even (odd) parity. The 
forms of the final expansions are given in Appendix ISl 



3.4. Spectral filter 

The geometry of the sphere makes grid points cluster naturally near the poles. In order 
to deal with the clustering, spectral methods use a filter to suppress high-wavenumber 



modes near the poles in the (p expansion ( 


Umscheid & Sankar Rao 


1971 


Fornberg & 


Merrill 


1997 


Bagchi & Balachandar 


2002 


Giacobello 


2005 


). From previous studies on 



the stability of swirling flow past a sphere, it is known that the k = ±1 modes {k is 



the azimuthal wave number) are the most unstable (Natarajan & Acrivos 1993). From 
the CFL condition, it can be deduced that the time step is determ ined by the A; = ±1 
modes if /3;jfe, "fijk, and Sijt decay faster than (lBagchi'12002 1. A filter that fulfills 



these conditions was devised by Bagchi fc Balachandar^ (2002p , in which the coefficients 
of the 4> expansion are multiplied by 



50(r,0,fc) = l-expK(fc)y^( 



fc)i 



(3.2) 



In (3.2 1, we define Y — rsm9, and ^(/c) and 0(fc) are functions subject to the boundary 
conditions gcf,{r,9i,k) = k~^ and g^{kr,9i,k) = 0.9. The exponential form of the filter 
ensures that its effects are limited to a small region near the poles of the sphere. Figure 
[TJi illustrates the behaviour of as a function of Y and fc. 

Aliasing arises because we are restricted to a finite range of wavenumbers ( Boyd|2001 1 . 
As a remedy, we adopt Orszag's 2/3 anti-aliasing ("padding") rule, which filters out 



waves with wavelengths twice and thrice the grid spacing. Orszag (1971&) showed that 
one obtains an alias-free computation on a grid with TV points by filtering out the high 
wavenumbers and retaining only 2A^/3 modes (BoydpOOT Canuto et a?.||1988 l. 

Spectral methods are sensitive to boundary conditions. Oscillations generated by the 
Gibbs phenomenon ( Boyd|[2001| contaminate the solution and grow unstably with time. 
In order to mitigate these instabilities, we multiply the coefficients of the r expansion by 
an expression of the form ( Voigt et a?.||1984 Don||1994 l 



cr (/, N.r) = exp — 



l-N, 



ln< 



(3.3) 



and the coefficients of the 9 expansion by a similar expression a{l,Ne), where ^ \l\ ^ Nr 
is the radial wave number, e = 2.2 x 10~^^ is the machine zero, 7 is the (integer) order of 
the filter, and N, is the central wavenumber of the filter. A small order (7 < 16) indicates 



8 

[htpb] 



C. Peralta, A. Melatos, M. Giacobello and A. Ooi 





0.15 0.20 



Figure 1. (a) Pole filter as a function of cylindrical radius Y and azimuthal wavenumber 
k, with 7 ^ fc ^ 24. The effects of the filter are greatest in a small region near the poles, 
whose cylindrical radius increases with increasing k. (b) Spectral filter a as a function of radial 
wavenumber I (or equivalently, latitudinal wavenumber j) and filter order 7, with 4^7^ 22. 
The filtering is weaker as 7 increases. The modes k = 0,±1 are not filtered, since they are 
important to get the correct stability characteristics. For all three wavenumbers, = 1 for all 
values of Y . 



strong filtering, while a high order (7 > 16) indicates gentle filtering. Figure [Tja shows 
the behaviour of the spectral filter (3.3) as a function of wavenumber, for 4^7^ 22. 



3.5. Initial and boundary conditions 

3.5.1. Initial conditions for v„ and Vs 

The velocity fields must be divergence-free initially, in order to satisfy the incom- 
pressibility constraint (2.3 1. The easiest choice is v„ = = 0. However, the superfluid 



velocity field is used to calculate the vorticity unit vector, cbs, which in turn appears in 



(2.6 1 and (2.7 1 and must remain well defined. Additionally, the HVBK equations describe 
a rotating superfiuid, implying a;^ 7^ in general. A simple initial condition that satisfies 
V • v„ = V • = and ^ is the Stokes solution ( ^Landau fc Lifshitz||1959 1 



0.1 



sin 9 Bj 



(3.4) 



This ansatz is an exact solution of the spherical Couette problem in the limit Re and 
a respectable approximation for Re < 10, where meridional circulation, which is absent 
from (3.4 1, carries only ^ 0.01 % of the total kinetic energy ( |Dumas|p9T| ). 



3.5.2. Boundary conditions for v„ 

The normal fluid satisfies a no-penetration condition, (v„)r = 0, at the inner and 
outer spheres. It also satisfies a no-slip condition; its tangential velocity equals that of 
the surface, like for a viscous, Navier-Stokes fiuid. The angular velocity vector which 
is tilted with respect to the z axis in the x-z plane by an angle can be written in 
spherical polar coordinates as 



$^2 — fl2 [(cos 9o cos 6 
— sin 6*0 sin (pe^ 



sin ^0 sin 9 cos <j))er — (cos 9o sin 9 — sin 6*0 cos 9 cos (j))eg (3.5) 
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while rJi remains fixed parallel to the z axis. The no- slip and no-penetration boundary 
conditions then reduce to 



v„(i?i,6',0) = RiVli X e^, 



(3.6) 
(3.7) 



3.5.3. Boundary conditions for 

The distribution of quantized vorticity in the superfluid component determines the 
boundary conditions for v,,. Quantized vortices in a cylindrical container are arranged in 
a rectilinear array parallel to the rotation axis if the rotation is slow [Re < 268; Barenghi 



& Jones (19881; Barenghi (1995)] or axisymmetric (Henderson et al. 1995 Henderson 
& Barenghi 2000 2004). Under these conditions, the numerical evolution is stable if 
the vortex lines are parallel to the curved wall (i.e. perfect sliding, x n = 0) and 
perpendicular to the end plates. 

In more general situations, e.g. noncylindrical containers, nonaxisymmetric flows, or 
fast rotation, there is no general agreement on what boundary conditions are suitable 
( Henderson fc Barenghi|2000 1 . This is especially true when the rectilinear vortex array is 
disrupted by the Donnelly-Glaberson instability to form an isotropic, turbulent vortex 
tangle (Section |2.2|. The radial component of the superfluid satisfies no penetration: 



(v,),(i?i,0,0) = O=(v,),(i?2,f?, 



(3.8) 



It is less clear how to treat the 9 and </> components. Numerical solutions of the HVBK 
equations in cylindrical Couette geometries are stable only if there is perfect sliding 



at the inner and outer surfaces (Henderson et al. 1995 Henderson & Barenghi 1995 
2004); numerical instabilities grow at rough surfaces ( Henderson et aZ.|1995 |. In spherical 



containers, however, the vortex lines are neither perpendicular to the walls nor parallel to 
the rotation axis everywhere. Previous attempts to solve the HVBK equations in spherical 



geometries foundered partly because of these issues Henderson & Barenghi (20041; R. 
Hollerbach 2004, private communication]. 



Khalatnikov ( 1965 ) suggested that vortex lines can either slide along, or pin to, the 



boundaries, or behave somewhere between these two extremes. If the boundary is not 



moving, the vortices terminate perpendicular to the surface (Khalatnikov 1965). The 



tangential velocity of a vortex line relative to a rough boundary moving with velocity 
u is given by ( |Khalatnikovl[T965l [Hills fc Robertsl[l977l [Henderson fc Barenghi|[2000| 

vl - u = ciWs X (n X Lbs) + C2n x w^, (3.9) 

where n is the unit normal to the surface, and ci and C2 are coefficients describing the 



relative ease of sliding. The form of (3.9) follows from calculating the energy dissipated 
as vortices slip along the surface ( Khalatnikov|196^ I . Equation (3.9) is difficult to include 
in HVBK theory, where each fluid element is threaded by many vortex lines, because 
is the velocity of a single vortex line; it cannot be calculated from v„ and v^, which are 
averaged over regions containing many vortex lines. Additionally, the slipping parameters 
Ci and C2 must be evaluated at each point on the surface, yet there is no experimental 
or theoretical study available in the literature on the precise form of these parameters. 
However, we can consider two simple limits of equation (3.9). For ci = C2 oo 
vortex lines slide freely along the surface and one requires 



the 



u)^ X n — 



(3.10) 



in order that remains flnite; that is, the vortex lines are oriented perpendicular to the 
surface. On the other hand, for Ci — C2 — 0, we have rough boundaries with = u. In 
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spherical Couette geometries, this implies no-slip, i.e. 



(3.11) 
(3.12) 



We find empirically that conditions (3.111 and (3.12) lead to stable numerical evolution 
in most scenarios studied in this paper. 

The existence of a vortex-free {ijJs — 0) region adjacent to the boundaries, whose 
thickness approaches the intervortex spacing, is theoretically predicted by minimizing 
the free energy of a vortex array in a container ( Hall] 1960 StaufFer fc Fetter]] 1968 Hills 
fc Roberts||1977| Henderson et a?.||1995 ). However, it has not been detected conclusively 
in experiments ( [Northby fc Donnelly 1970 Mathieu et a?.|1980 |. It is unclear how to treat 
this boundary layer numerically within HVBK theory, which assumes a high density of 
vortices, so we do not consider it further in this paper. 

3.5.4. Accelerating spheres 

The angular velocities of the outer and inner spheres, VI2 and rJi, can vary with time, 
either in a prescribed way or in reaction to the torque exerted by the fluid. 

One scenario considered in this paper is free precession of the outer sphere. This 
situation is relevant to astrophysical systems like neutron stars ( Jones fc Andersson|2002 



Shaham|1977 Link 2003 Sedrakian|2005 ) and to laboratory systems like superfluid-fiUed 
gyroscopes (Reppy 19651. Let the outer sphere be biaxial, with symmetry axis 63 and 



constant total angular momentum J, and resolve the angular velocity into components 
( |Shaham|[T986t 



(3.13) 



Jii Jj^ 

parallel and perpendicular to the symmetry axis, with J = J|| + Jj^, where I\\ and /_l are 
the associated moments of inertia. The precession frequency ilp is then defined by 



1 



1 



(3.14) 



and the velocity of any point on the surface of the outer sphere in the inertial (lab) frame 
is given by 

— — R2^' X Br — R2^pe^ X Br, (3.15) 

where O! — il'3/J is the inertial- frame precession frequency. The back-reaction of the 
fluid on the container needs to be included when solving the HVBK equations self- 
consistently. The viscous torque accelerates (decelerates) the container. To this must be 
added any external torques iVext- Examples of TVext are the electromagnetic torque on the 



crust of a neutron star (Ostriker & Gunn 1969 Melatos 1997 Spitkovsky 20041 or the 



friction between a rotating container and its spindle in laboratory experiments (Tsakadze 
fc Tsakadze||1972 19801. In this situation, r2i and ^2 evolve according to 

1,2 



dn 



dt 



(3.16) 



where lij is the moment-of-inertia tensor, and N^^^ is the instantaneous viscous torque 
exerted by the normal fluid on the shell ( Landau fc Lifshi"tz||1959 1, 



Re 



dO r sin ( 



dVnd 



Vn<P 

r 



(3.17) 
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Equation (3.16 1 is solved explicitly at each time step using a third-order Adams-Bashforth 
algorithm to get ^i{t + At), after advancing the flow using ni{t) (see Section 3.2 and 
Appendix C.l ). 



4. Validation 



To the best of the authors knowledge, the problem of superfluid SCF has never been 
solved before for Re ^ 1, save for an inconclusive pioneering attempt by R. HoUerbach 



[private communication, 2004; see also Henderson & Barenghi (2004l], who encountered 
numerical instabilities when implementing the cylindrical boundary conditions used by 
IHenderson et al. ( 1995 1 . Consequently, we cannot verify our code directly against previous 



superfluid SCF results, and we are forced into a different validation strategy: in the limit 
T Tc, the superfluid component vanishes, equation (2.11 reduces to the classical Navier- 



Stokes equation, and we validate our numerical scheme against the wealth of numerical 
and experimental studies available for viscous SCF. 

Our three-dimensional pseudospectral HVBK solver reduces to a classical Navier- 
Stokes solver if all the coupling terms [HVBK friction, HVBK tension, V(v^5) and 
z^sVws] are removed, and the arrays are disabled. We validate our solver for three 
types of SCF in this regime: (i) inner sphere rotating, outer sphere stationary (fii ^ 0, 
^2 = 0); (ii) inner sphere stationary, outer sphere rotating (fii = 0, 0); and (iii) 

both spheres rotating (fii ^ 0, 7^ 0). For the parameter range explored in this Section, 
viz. 50 < Rei,Re2 < 1200 and 0.18 < (5 0.5, a grid with {Nr,Ng,N^) = (81,200,4) 
and At = 10~^ is sufficient to fully resolve the ffow. A ffow is regarded as fully resolved 
if the spectral mode amplitudes decrease quasi-monotonically with polynomial index. 

The meridional streamlines drawn in the flgures below correspond to the final steady 
state. A steady state is deemed to have been reached when the difference on the viscous 



torques between the inner and outer spheres satisfies 
expressed in units of pR^il^ or pR\Q. 



^cxt - Nit 
2, unless otherwise indicated. 



^10 Torques are 



4.1. Inner sphere rotating, outer sphere stationary 



Following Marcus & Tuckerman (1987a I, we focus on a moderately sized gap 5 — 0.18 



with 0,1 ^ Q and 0,2 = 0. We look for time-dependent transitions between axisymmetric 
steady states characterized by zero, one, or two Taylor vortices on either side of the 
equatorial plane, and refer to them as 0-, 1-, and 2-vortex states respectively. We obtain 



these basic flow states, together with an intermediate "pinched" vortex state (Bonnet &: 
Alziary de Roquefort||1976 ), as illustrated in Figures 4, 7 and 9 of |Marcus fc Tuckerman 
( |1987aD . 

Figure [2] plots the steady-state viscous inner torque (3.17) as a function of Reynolds 
number for 50 ^ Rei ^ 600, normalized to the torque exerted by Stokes flow, A^stokcs = 
167r/[i?ei(l — -Rf/i?!)] (Marcus & Tuckerman 1987a I. The square symbols record the 



values taken from Figure 1 of Marcus & Tuckerman ( 1987a'), while the circles are output 



from our numerical code. Each poin t is obtained by starting from an initially stationary 



fluid, not the Stokes solution (3.4), and evolving in time until a steady state \N- 



(2) 



^10 is reached. The points obtained from our numerical simulations and the 



published results agree to three signiflcant digits. 

Once the basic 0-vortex state is attained for Rei = 600, transitions to pinched, 1- 
vortex, and 2-vortex states can be induced by impulsively changing Rei (by reducing the 



viscosity) to 650, 700, and 900 respectively. Below we check against Marcus & Tuckerman 
( 1987a|& I whether our solver follows these transitions faithfully. 



12 

[htpb] 



C. Peralta, A. Melatos, M. Giacobello and A. Ooi 



1.06 



1.04 



1.02 




600 



Re 



Figure 2. Viscous inner torque on the inner sphere, in units of the Stokes torque 
A^stokc. = 167r/[7?ei(l - Rf/R^)], versus Reyolds number Re. T he circles correspond to outpu t 
from our numerical solver. The squares denote data taken from [Marcus fc Tuckerman| ( |1987ffll ). 



4.1.1. 



1 transition 



We simulate the — > 1 transition by starting with a 0-vortex equilibrium at Re = 
650 and then abruptly (over one time step) reducing the viscosity to give Re = 700, 



where the equilibrum becomes unstable (Marcus & Tuckerman 1987&I. We obtain the 



intermediate states displayed in Figure 10 in Marcus & Tuckerman (1987&I. At the start 



of the sequence, the streamlines are not symmetric about the equator; the boundary 
between the counterrotating vortices at t = 230 is displaced south of the equator. Then, 
at t = 315, two wedges start to form in the northen hemisphere, at ~ 82 deg, and generate 
a growing vortex ai t — 320, which evolves into a fully developed Taylor vortex in the 
northen hemisphere at i = 400. The transitions occur at about the same time as Figures 



lOa-lOb in Marcus & Tuckerman 



lOc-lOf in Marcus & Tuckerman 



( 1987& ), and about 68 units of time later than in Figures 
(1987&I. Figure [3] shows how the torque evolves during 



the time interval covered by the 0^1 transition. The transition, marked by a jump 
in the torque, occurs later than in the numerical experiments of [Marcus fc Tuckerman] 
(1987&I because it is sensitive to the exact form of the initial perturbation, which in 



turn depends on roundoff error, aliasing, and resolution (Marcus fc Tuckermanl 1987&I 



The sensitivity is exacerbated by the north-south asymmetry of the 0^1 transition (cf. 
— > 2 below). However, the shapes of the solid and dashed curves in Figure [S] (especially 
the growth rate) agree to three significant digits if we slide them on top of each other. 



4.1.2. — + 2 transition 



According to Marcus fc Tuckerman (1987&I, the 0^2 transition can be produced by 



starting with a 0-vortex equilibrium (50 < ^ei < 651) and impulsively increasing Rei 
above 775, where the vortex equilibrium is unstable. We therefore start with the Stokes 
solution for Rei = 650 and suddenly increase Rei to 800. We obtain the 0^2 transitions 





I 

200 




400 



60U 



SOU 



Time 

Figure 3. Torque during tlie 0^1 transition, (a) Torque on the inner sphere, A'^i, as a function 
of time, (b) Difference between the inner and outer torques as a function of time. The dashed 
curves correspond to numerical results from Marcus & Tuckerman ( 19876 I, while the solid curves 
are generated by our numerical code. 



of the meridional flow as in Figure 4 by Marcus & Tuckerman (1987&). In Figure |4] we 
plot the torque on the inner sphere and the difference between the inner and outer torques 



as functions of time (solid curve), together with data taken from Figure 4 of Marcus & 
Tuckerman (1987&J (dashed curve), showing an agreement to three significant digits, 
after sliding the curves together. In this case the transition is symmetric with respect 
to the equator and occurs more quickly. The — s- 2 transition is always symmetrical 
about the equator, as compared to the — s- 1 transition presented in Section 4.1.1 In a 



bifurcation diagram (showing the relation between torque and critical Rei), the 0- vortex 
and 2-vortex flows lie on the same critical branch, while the 1-vortex state lies on a 
different, non intersecting branch ( Marcus fc Tuckerman|[T987& I. 



4.2. Outer sphere rotating, inner sphere stationary 

We now allow the outer sphere to rotate (^2 ^ 0), while holding the inner sphere fixed 
(fii = 0). We reproduce the meridional streamlines for Re2 — 100, 500, 1000, and 2000, 
with 6 — 0.5, obtained by Dennis fc Singh] ( 1978| ) (Figures 3, 4, and 5 of their paper), who 
used a numerical method in which the flow variables are expressed as a truncated series 
of n orthogonal Gegenbauer functions with variable coefficients, reducing the Navier- 
Stokes equation to a set of ordinary differential equations. At Re2 = 500, the agreement 
is good, with the streamlines showing the same distribution of vortices in the northern 
hemisphere: one primary circulation cell, slightly elongated in the direction of the rotation 
axis, with its center located ~ 40 deg over the equator, and a small recirculation zone 



14 

[htpb] 



C. Peralta, A. Melatos, M. Giacobello and A. Ooi 



(a) 



G.lJ 




I , , I 

Ell] ■,■[) ii\ !HI 1(1(1 111) 



(b) 




_L 



60 
'I 'iirif 



j: 



"JO so W) JOO 110 120 



Figure 4. Torque during the 0^2 transition, (a) Torque on the inner sphere, A^'i, as a function 
of time, (b) Difference between the inner and outer torque as a functio n of time. The dashed 
curves correspond to numerical results from Marcus & Tuckerman ( 1987a I, while the solid curves 
are generated by our numerical code. 



near the equator. For Re2 — 1000, the primary circulation cell is more elongated, with 
most of the circulation lying in a cylindrical sheath of radius approximately equal to 
as predicted by the Taylor-Proudman theorem ( Pearson]|1967 1. For Re2 — 2000, the 
agreement is not as good. Dennis fc Singh, ( ,1978| were unable to obtain a well defined 
flow pattern for Re2 — 2000, having been limited by computational resources to only 
eight Gegenbauer polynomials, which is not sufficient to follow small vortex structures 
developing near the equator. The higher-resolution results in our simulations suggest 



that the small vortices observed by Dennis & Singh (19781 near the equator are probably 
low-resolution artifacts. 

The steady-state dimensionless torque calculated by various authors (including the 
present work) is presented in Table [l] together with bibliographic information. 

4.3. Inner and outer spheres rotating 
The next step in the verification program is to consider the rotation of both spheres. We 



follow Pearson (19671 and Munson & Joseph (19711, who studied general axisymmetric 



flows between rotating spheres, solving the Navier-Stokes equation in terms of stream 



functions in a meridional plane. Pearson (1967) used a numerical scheme based on finite 



differences (with typical resolution of 40 x 20 mesh points) , whereas Munson & Joseph 
(19711 used expansions in Legendre polynomials (typically using up to 7 terms). 

Suppose the spheres counterrotate, with fii = ~^2- We obtatin the meridional stream- 



lines and angular velocity profiles as shown in Figures 4-5 of Munson & Joseph ( 1971 1 
for Re2 = 100 and 500 respectively, with 6 = 0.5. For Re2 = 100, the agreement is good, 
with the 0-vortex state rotating clockwise in the northern hemisphere, counterclockwise 
in the southern hemisphere, and centered at « 45 deg. The angular velocity contours are 
nearly concentric shells with values decreasing from Q.2 at the outer shell to —^2 at the 
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Re2 N2 [pRln 



Reference 



100 



500 



0.041745 
0.041750 
0.042450 
0.041888 
0.04160 



Dumas 



Present st udy 

( 1991|); Dumas fc Le onard | (|1994 1 



Dennis & yuartapcllc ( 198^ 
Dennis Mmsh 



Munson 



&«ingh(1978J 
& JosepliTraTl 



0.011979 
0.011985 
0.012282 
0.011980 



P resent study 

Dumas' fl991') ; Dumas & Leonard 



Ucnnis fc Cjuartap ellc (1984| 
' — |Deiirn7fcSiiiih| |I5f§] r 



(19941 



1000 7.2203 X 10 
7.2237 X 10 
7.7074 X 10 
7.2382 X 10 



P resent study 

Dumas' ('1991'); Dumas & Leonard 
Dennis & Quartap clle (^1984 ) 
' iDennis & Singh) (pfSt—^ 



(19941 



2000 4.4483 x 10 
4.4478 X 10 



Present stud y 



Dennis & Singh ( 1978 1 



Table 1. Comparison of numerical values obtained by various authors for the torque on the 
outer sphere, N2, when the outer sphere is rotating and the inner sphere is stationary. 



inner shell. At Re2 — 500, an additional counterclockwise vortex develops in the polar 
region, near the inner sphere, because the influence of the inner sphere strengthens as 
the viscosity decreases. The locations of this vortex and the main circulation cell (with 
its center at « 30 deg) agree with the results of Munson & Joseph ( 1971 1. The angu- 
lar velocity profiles show a similar pattern, forming a cylindrical sheath parallel to the 
rotation axes. 

Now suppose that the spheres counterrotate, but with fii = —2^2- We do numerical 
simulations for Rt — U1R2/1' = 100 and Re = 500 respectively. We get good agreement 



with the simulation results of IMunson & Joseph (19711 presented in Figures 6 and 7 of 



their paper. The faster rotation of the inner sphere produces an additional circulation 
cell near its surface, both for Re = 100 and Re = 500. The center of the secondary cell is 
slightly displaced towards the equator in the latter case. The angular velocity contours 
tend to form a cylindrical sheath as the Reynolds number increases ( Proudman||1956 1. 

Figure |5] plots the dimensionless torque as a function of the Reynolds number Rei 
or Re2 (the definition used in each case is indicated in the plots). When the inner and 
outer spheres rotate in opposite directions, with fii = — 2f22, the inner torque is shown in 
Figure [5^; for fii — —U2, the inner torque is shown in Figure [5|d. When a steady state is 
reached, the difference between the inner and outer torques approaches ~ 10~^pi?2^f 2- 



The solid curve and asterisks represent the data obtained by Munson & Joseph ( 1971 1 



and from our numerical simulations respectively. The results coincide to three significant 
digits for all the Reynolds numbers considered when fli = —2172. However, the results 
diverge for Re ^ 500 when fii = —U2- This is a consequence of the low resolution in the 



expansions used by Munson & Joseph (1971), who claimed to be unable to reproduce 
small structures near the equator, of typical size ~ 0.3i5, when comparing with the study 
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0.9 



1000 
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Figure 5. (a) Dimensionless torque on the inner sphere, Cm = 7 ReN\ / 1Q-k pB^Q.\, as a function 
of Reynolds number, Re. The spheres rotate in opposite directions, with Q,i — —2Q2- (b) Dimen- 
sionless torque on the inner sphere. Cm = 2ReNi/3TvpR2Q,2 as a function of Reynolds number. 
Re. The spheres rotate in opposite direct ions, with f^i = —Q2- The solid curve corresponds to 
data taken from Munson & Joseph (19711. The asterisks are output from the present study. 



by Pearson (1967). Munson & Joseph (1971) used a maximum of 7 modes in r and 9, 
whereas we use A^^ = 81 and Ng = 200 and can therefore resolve vortical structures as 
small as 10~'^6. Note that the torque is dominated by surface regions where the shear 
stresses are stronger, e.g. where vortices cluster. 



5. Unsteady, superfluid SCF 

In this section, we investigate the unsteady behaviour of SCF in classical (Navier- 
Stokes) fluids and superfluids in two dimensions, by performing a set of axisymmetric 
numerical experiments {N^ — 4) with rotational shear in the range 0.1 ^ Ail ^ 0.3, 
in medium and large gaps (0.2 ^ 5 ^ 0.5). For HV mutual friction, we use B = 1.35 

Donnelly||2005j 



and B' = 0.38, the He II values at T = 1.45 K (Barenghi et al. 1983 



Donnelly & Barenghi 1998[ ). [fjFor GM mutual friction, the parameter A' = 5.8 x 10 



t We consider adiabatic walls and divergence- free Vs and v„. Although one expects the tem- 
perature to rise continually in this scenario due to dissipation, we ignore the influence of dis- 
sipation inhomogeneities in the superfluid flow. This is equivalent to assuming psV^ <C PnC^, 
where v — Vn + Vs and c is the second sound speed. In all the simulations presented in this 
paper we have O.Olv^ < 0.99c^ We can calculate the rate of change of the internal energy of a 
unit mass of fluid due to viscous heating from E = Rea^^A/2pn, where A is the total time of 
the simulation (A ~ 10^) and cr,.0 is the viscous stress tensor 
parameters used in the simulations, we have 10-* < S < 10" 



Landau fc Lifshitz| ( |1959[ ). For the 
'', which is safe to ignore . 
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Figure 6. Snapshots at t = 214 of meridional streamhnes for the normal (left) and super- 
fluid (right) components in superfluid SCF, with Re = 10*, 5 — 0.5, and Af2 — 0.3. Spectral 
resolution: (Nr,No,N^) = (150,400,4). Filter parameters: (7r,7e) = (8,6). 



(with Xi/X2 = 0.3) at the same temperature can be calculated from a fitting formula 
derived by Dimotakis (19721, which is consistent with previously published experimental 
values ( Vinen|1957 



Stable long-term evolution is difficult to achieve for this value of A', 
5.8 X 10~^ instead. We compile a preliminary gallery of vortex states, in 



so we take A' 

the same spirit as for classical SCF and the validation experiments in Section |4] ( ,Marcus 
fc Tuckerman||1987a|& Dumas ||T991 Junk fc Egbers|[2000 ); a complete classification lies 



outside the scope of this paper. The torque is observed to oscillate quasiperiodically, yet 
persistently, accompanied by oscillations in the vortical structure of the fiow (Sections 



5.1 5.2). Resolution and filtering issues are discussed in Section 5.3 The role played by 



the inviscid superfiuid, and the effect of varying the strength and form of the mutual 
friction force, are studied in Section |5.4| . It is observed that the superfluid tends to 
destabilize the flow and increases the torque. Boundary conditions are varied in Section 
1531 

5.1. Meridional streamlines 

Figure |6] depicts the meridional streamlines of the normal (left) and superfluid (right) 
components in superfluid SCF, for the special case Re = 10**, S = 0.5, and Ail — 0.3. In 
the equatorial zone (60° <0< 120°), we observe two large circulation cells adjacent to 
the inner boundary. Each large cell contains twin cores circulating in the same sense (and 
therefore tending to repel). Between the large cells and the outer boundary exist two small 



18 



C. Peralta, A. Melatos, M. Giacobello and A. Ooi 



vortices, occupying « 20 % of the volume of the large cells. The flow in each hemisphere 
is symmetric about the equatorial plane. Away from the equator (30° < |0 — 90°| < 90°), 
we observe a large cell (width w 30°) at mid latitudes, three vortices in the normal 
component, and one small vortex in the superfluid component. 

The flow pattern described in the previous paragraph is characteristic of moderately 
high Reynolds numbers (Re > 10*). The HV mutual friction couples normal and su- 
perfluid components strongly, so that their meridional streamlines are similar. At lower 
Reynolds numbers (-Re < 10"^), the streamlines of the two components differ markedly. 
The normal component behaves like a viscous, Navier-Stokes fluid at low Re, with a 
small number (< 3) of large circulation cells on each side of the equatorial plane. The 
superfluid is influenced less by the normal fluid, due to the stiffness provided by the vor- 
tex tension force ( Henderson fc Barenghi|1995 Swanson|1998 ). Streamlines of Vs develop 



multiple eddies and counter-eddies. When GM mutual friction operates, the normal and 
superfluid components behave similarly, both at low and high Reynolds numbers, but the 
superfluid displays a richer variety of circulation cells, while the normal component be- 
haves like an uncoupled Navier-Stokes fluid. The different effects of HV and GM mutual 
friction are investigated in Section [5.4[ 

Figures [t] and [s] show the meridional streamlines ai t — 20 and Re = 3x10* for 
the normal and superfluid components, with GM friction and zero tension force. The 
rotational shear 0.1 ^ Ail ^ 0.3 increases from left to right; the dimensionless gap width 
0.3 ^ (5 ^ 0.5 increases from top to bottom. The flow is approaching, but has not reached, 
a steady state at this time, with \N^'^l — N.^^^\ < 10~'^. The number for equatorial and 
polar circulation cells remains approximately constant as 6 increases, although the cells 
become progressively less "stacked". By contrast, the flow becomes less turbulent, and 
the number of cells decreases. Additionally, the vortices are more stacked with decreasing 
S. 

Figures [9| and [TO] show meridional streamlines at i = 14 for the normal and superfluid 
components as a function of Reynolds number, for fixed S and Af2, HV friction, and 
nonzero tension. The streamlines of the normal fluid show the 0-vortex state at Re — 100 
(Figure |9|i) . A secondary vortex develops near the outer shell at Re — 300 (Figure ^p) , 
whose size increases with Re, elongating in the meridional direction. Two additional 
vortices form in the equatorial region at Re — 10^ (Figures [9|;-f). The streamlines of the 
superfluid component closely resemble the normal fluid component for Re > 3 x 10'^ (see 



Figures lOi-f) 



5.2. Unsteady torque 

The torque exerted by the normal fluid component on the inner and outer spheres, 
calculated using (3.171, is plotted versus time in Figures 111 and[TT[3. It oscillates, with 



peak-to-peak amplitude ~ lO^'^ for i ^ 30 and ~ 10^^ for t ^ 30. These oscillations, with 
period « 27r/rii, persist as long as the differential rotation is maintained, up to t = 214 
in our longest simulation. They are observed at all the Reynolds numbers considered in 
this paper (1 x 10^ ^ Re ^ 3 x 10*). The oscillation amplitude is greater for HV friction; 
oscillations are still observed for GM friction, but with peak-to-peak amplitude ^ 10~^. 
In other words, superfluid SCF is intrinsically unsteady and indeed quasiperiodic. 

5.3. Spectral resolution and filter 

In a well resolved simulation, the Chebyshev and Fourier mode amplitudes decrease 
monotonically (on average) with polynomial index. We prefer to maintain an amplitude 



ratio of at least 10 between the strongest and weakest modes. Giacobello (20051 found 



empirically that this is sufficient to fully resolve vortical structures in unsteady, three- 
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Figure 7. Snapshots at f = 20 of meridional streamlines for the normal fluid component in 
superfluid SCF, with _Re = 3 x fO"*, 5 = 0.3, 0.4, 0.5 (from top to bottom), and AO. = 0.1, 0.2, 
0.3 (from left to right). The friction force if of GM form, with zero tension (T = 0). Spectral 
resolution: {Nr,No,N^) = (120,250,4). Fiher parameters: (7r,7e) = (8,8). 



dimensional transients excited by the flow past a stationary and rotating sphere in a 
classical, viscous Navier-Stokes fluid. In this section, which is devoted to axisymmetric 
flows, we are interested in the Chebyshev (r) amplitudes w^r,sr; "^raesSi ^^'^ ''^ruf) s4>i ^"^^ 
Fourier {9) amplitudes w^^.sr' "ne sS' ^^'^ ''^ncp s4>- These coefficients do not correspond to 
Pijk, li]k, and 6ijk in equations (B2)-(B4|; v]^^ ..^, ^ne^se^ ^^'^ K^,s,p ^re calculated by 



transforming from coordinate to wavenumber space in the r direction, and 

and are calculated by transforming from coordinate to wavenumber space in the 

6 direction. 

We start by comparing the mode amplitudes for a poorly resolved and a well resolved 
simulation. Figure 12 shows an example of a poorly resolved simulation, with Re = 178, 
6 = 0.5, An = 0.1, and {Nr,Ne,N^) = (150,400,4). The spectral fllter is switched off. 
The mode amplitudes of the normal velocity components decrease from ~ 10^^ (i = 1) to 
~ 10~^ {i = 150), which superficially suggests that the fiow is well resolved. However, the 
mode amplitudes of the superfluid velocity components are roughly constant (~ 10~^) 
across all Chebyshev (1 ^ i ^ 150) and Fourier (1 ^ j ^ 400) orders, indicating that 
the run is poorly resolved. Figure [13] shows a well resolved simulation for the same 
parameters. The improvement is achieved by switching on the spectral filters, and the 




Figure 8. Snapshots at t = 20 of meridional streamlines for the superfluid component in 
superfluid SCF, with Re = 3 x lO", S = 0.3, 0.4, 0.5 (from top to bottom), and AO. = 0.1, 0.2, 
0.3 (from left to right). The friction force if of GM form, with zero tension (T = 0). Spectral 
resolution: {Nr,Ng,N^) = (120,250,4). Fiher parameters: (7r,7e) = (8,8). 



extent of the improvement depends on the strength of the filters, with 7,, = 8 and jg — 6 
in Figure [13] The Chebyshev amplitudes of v„ decrease by 8 orders of magnitude over 
the index range 1 ^ i ^ 25. The Chebyshev amplitudes of Vg decrease more gradually, 
by '--^ 8 orders of magnitude over the range 1 ^ i ^ 100. The Fourier amplitudes of v„ 
and Vs decay similarly (Figures [Tsji-d) , with only Ng sa 20 required to resolve the flow 
properly. For a weaker filter, with 7,, = 79 = 12, the mode amplitudes are unchanged to 
within ~ 10 %, but Nr « 20 and Ng « 100 Chebyshev and Fourier modes are required. 

What happens in general when the exponential filter is either switched off, as in Figure 
12 or maintained at a weak level (7,., 75 > 20)? For Re > 10^, we find that the evolution 
IS stable for a short time {t < 10), after which Vji and become unphysically large and 
the numerical simulation breaks down. For Re ^ 10^, the evolution is stable for longer, 
provided that perfect slip boundary conditions are apphed to v^. Indeed, generally speak- 
ing, Vg evolves less stably for no-slip boundary conditions, which promote the generation 
of superfluid vorticity. Nevertheless, for a range of SCF parameters, we observe that 
v„ for a flltered HVBK superfluid agrees well with v for an unflltered Navier-Stokes 
fluid given identical boundary conditions (see Section 5.4 1, engendering confidence that 
filtering does not cause artificial behaviour. 
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Figure 9. Snapshots at f = 14 of meridional streamlines for the normal fluid component in 
superfluid SCF, with 5 = 0.5, Af2 = 0.3, and Re increasing from top left to bottom right, (a) 
Re = 100, (b) Re = 300, (c) Re = 1000, (d) Re = 3000, (e) Re = 10*, and (f) i?e = 3 x 10"*. 
The friction force is of HV form; the coefficient of the tension force is Rcs ~ 10~^. Spectral 
resolution: {Nr,Ne,N^) = (120,250,4). Fiher parameters: (7r,7e) = (8,8). 
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Figure 10. Snapshots at t = 14 of meridional streamlines for the superfluid component in 
superfluid SCF, with S = 0.5, AQ = 0.3, and Re increasing from top to left to bottom right, (a) 
Re = 100, (b) Re = 300, (c) Re = 1000, (d) Re = 3000, (e) Re = lO", and (f) Re = 3 x 10*. 
The friction force is of HV form; the coefficient of the tension force is Rbb = 10~^. Spectral 
resolution: (Nr,N0,N,i,) = (120,250,4). Fiher parameters: (7r,7e) = (8,8). 



5.4. Effect of superfluid component 

Laboratory experiments on the acceleration and deceleration of He II in spherical vessels 
reveal a variety of unsteady behaviour, e.g. sudden jumps and quasiperiodic oscillations 
in angular velocity (Tsakadze & Tsakadze 19801. It is not known what aspects of this 
unsteady behaviour are caused by the nonlinear hydrodynamics of the viscous normal 
component of He II, or by the build-up of vorticity in the inviscid superfluid component. 
We now explore this question. 

In order to isolate how the superfluid component influences the global dynamics of su- 
perfluid SCF, we compare three particular cases: (i) a one-component, viscous, Navier- 
Stokes fluid, with Reynolds number Re — 10''; (ii) a one-component, inviscid, HVBK 
superfluid, with F = T = in (2.1) and (2.2); and (iii) a two-component, HVBK su- 
perfluid, whose normal component {Re — ICF) is coupled to the superfluid component 
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Figure 11. Viscous torque exerted on the (a) inner and (b) outer spheres as a function of 
time in superfluid SCF, with HV mutual friction, S = 0.5, Af2 = 0.3, and Re = 10^. Spectral 
resolution and filter parameters are the same as in Figure [6] Initially, we set v,i = Vs to the 
Stokes solution (3.4 1. 



through HV pall fc Vineii||1956a|&| ) and GM ( |Gorter fc Mellink|p49| |Donnelly|[2005 1 
mutual friction (F ^ 0, T ^ 0). To make the comparison, we fix (5 = 0.5 and Afi — 0.3. 

Let us begin with case (i): a viscous, Navier-Stokes fluid. Meridional streamlines, 
obtained by integrating the in-plane components of the velocity field in the x — Q plane, 
are displayed in Figure 14 1, at time t = 6. The flow is complicated, featuring three 
primary circulation cells: one near the equator, and two near the poles (each about half 
the diameter of the equatorial cell). Two secondary, flattened vortices reside near the 
outer boundary, one just above the equator and the other centred at r « 0.9, 9 ~ 30°. 
These structures develop from a low-J?e Stokes flow at t = 0. 

Let us repeat the experiment with an HVBK superfluid in which the normal and 
superfluid components are completely uncoupled, i.e. the mutual friction and tension are 
switched off (F = T = 0). This is case (ii). Naturally, the normal component evolves 
exactly like the classical Navier-Stokes fluid in Figure [141 fo'" the superfluid, its 
meridional streamlines at time t = 6 are shown in Figure [14|3. On large scales, the flow 
pattern resembles Figure 14 i, with the same number (three primary plus two secondary) 
of circulation cells. However, the cells have different shapes and diameters: the secondary 
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Figure 12. Snapshot of spectral mode amplitudes at t = 0.3, as a function of the Chebyshev 
polynomial index i for (a) the normal velocity resolutes (v„)r, (v„)a, and (v„)0, and (b) the 
superfluid velocity resolutes (va)^, (vs)e, and (vs)</,; and as a function of the sine (cosine) polyno- 
mial index j for (c) the normal velocity resolutes (v„)r, (v„)9, and (vn)^, and (d) the superfluid 
velocity resolutes (ys)r, iys)B, and (vs)^. The grid resolution is {Nr, Nq, N^) = (150, 400, 4) and 
the exponential filters are switched off. Simulation parameters: HV mutual friction, Re = 178, 
S = 0.5, and AO. = 0.1. 



cells are smaller and less flattened than in Figure [T4k, and the volume of the largest 



(equatorial) primary cell is three times greater than in Figure 14 1. 

Note that equations (1.16) and (1-17) evolve independently when the coupling forces 
(T and F) are switched off. The superfluid equation of motion (1.17) reduces to the 
Euler equation for an ideal (inviscid) fluid, so only one boundary condition is required: 
no penetration, {vs)r = 0. The same is true for HV and GM mutual friction: the two 
forms of the friction force imply two different orders of the system of equations and hence 
two different sets of boundary conditions. However, there are three reasons why we do 
not treat HV and GM friction differently in our simulations: 

(a) The correct boundary conditions for the superfluid are unknown. What ultimately 
determines the boundary conditions for Vj is the distribution of quantized vorticity in the 
superfluid component. In cylindrical containers, it is natural to assume that the quantized 
vortices are arranged in a regular array parallel to the rotation axis, if the rotation is slow. 
Under these conditions, the numerical evolution is stable if the vortex lines are parallel to 
the curved wall. In spherical containers, the vortex lines are neither perpendicular to the 
walls nor parallel to the rotation axis everywhere. In this paper, we test what boundary 
conditions give the most stable evolution. Often, no slip in the superfluid component 
works best, even if it is not strictly mathematically correct for the GM force. Physically, 
we justify this by assuming that vortices pin to the boundaries, whereupon they move at 
the same speed. In other words, no slip corresponds to pinning at the boundaries; since 
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Figure 13. Snapshot of mode amplitudes ai t = 0.3 as a function of the Chebyshev polynomial 
index i for (a) (v„)r, (v„)6i, and (vn)^, and (b) (vs)r, (vs)e, and (vs)^. Snapshot of mode am- 
plitudes a.t t = 0.3 as a function of the sine (cosine) polynomial index j for (c) (v„)r, (v,i)e, and 
(v„)0, and (d) (vs)r, (v,)e, and (v,)^. The grid resolution is {Nr,Ne,N4,) = (150,400,4). The 
exponential filter has (71, 7e) ~ (8,6). Simulation parameters; HV mutual friction. Re = 178, 
S = 0.5, and AO. = 0.1. 



we ignore the presence of the vortices in the fluid interior (T = F = at i?i < r < R2) 
but not at the boundaries. 

(&) When we attempt to solve the the HVBK equations numerically with T = in 
the HV friction force, or with GM friction, using only one boundary condition on the 
superfluid component (no penetration), the results do not differ significantly from the 
results presented in this paper. For 6 = 0.5, Ail — 0.3, at t = 4, the outer torque differs 
by 10~^ between the two approaches. One cannot see any significant difference in the 
streamlines of both fluids, whether one uses the strictly mathematically correct set of 
boundary conditions for the superfluid (no penetration only) or the simulations presented 
in this paper (no slip in the superfluid component). 

(c) The tension of the vortex lines, controlled by the stiffness parameter Vs, can dis- 
rupt the flow and destabilize its evolution. A high value of Vg {Rss < lO**) stiffens the 
vortex array against the drag of the normal fluid. Previous attempts to solve the HVBK 
equations in spherical geometries foundered partly because of these issues. By enforcing 
no-slip boundary conditions on the superfluid, one ensures that the superfluid corotates 
with the normal fluid at the boundaries, reducing the friction force in that region. As 



Henderson et al. ( 1995 1 suggested (page 333, second last paragraph), the appearance of a 



nonlinear boundary layer may not be completely resolved by our code, or indeed by any 
code to date. Note also that Vs is small 10~^), so the influence of the tension term 
in the simulations presented in this paper is not strong. We plan to study the effects of 
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the tension force in a future paper, when we are able to obtain a better understanding 
of the nonhnear boundary layer. 

When the mutual friction and tension forces are switched on, as in case (iii), the vortical 
structures near the poles change. Consider first what happens when F is of HV form. For 
the normal fluid component, displayed in Figure 14:, the flattened vortex near the pole 
widens radially to about twice the size of the same vortex with F = 0, while its latitudinal 
extent remains unchanged. The vortex near the equator halves in size, compared to its 
counterpart in Figu re |14^ . The only appreciable changes in the superfluid flow pattern, 
displayed in Figure |14|i, are the following: (i) the small vortex near the pole widens 
radially by ^ 10 %; and (ii) the larger circulation cell near the pole comes to resemble 
the same cell in the normal component more closely. 

By contrast, when F is of GM form, the streamlines of the normal fluid are very 
similar to a classical viscous fluid at the same Reynolds number, as we can appreciate by 
comparing Figures [T4| l and |14[ 3. The superfluid component closely resembles an uncoupled 
superfluid (c.f. Figures 14 3 and 14'). The two components are lo osely coupled because 
GM friction is weaker than HV friction [|Fhv/Fgm| ^ 10^; see Peralta et al. (2005 
2006&I I]. 



The evolution of the torque on the inner and outer spheres is plotted in Figure [15] for 
cases (i) and (iii). [The torque is zero in case (ii), where the fluid is completely inviscid.] 
The solid curve corresponds to a viscous, Navier-Stokes fluid, while the dashed and dotted 
curves correspond to an HVBK superfluid with HV and GM mutual friction respectively. 
Note that, in an HVBK superfluid, the torque is still exerted by viscous stresses in the 
normal component. However, its magnitude is modified by the presence of the superfluid 
component, because the mutual friction modifies v„ and hence d{'Vn)i/dxj. 

The torque exerted by the normal component increases roughly thrice relative to case 
(i) when the superfiuid component is included with HV mutual friction. By contrast, the 
torques exerted by a Navier-Stokes fiuid and a HVBK superfiuid with GM mutual friction 
differ by 6 %, which is barely distinguishable on the scale of Figure 15 To understand 
this effect quantitatively, consider the streamline snapshots of the Navier-Stokes fiuid, 
viscous HVBK component, and inviscid HVBK component at i = 20, shown in Figures 



16 1- 16 There are four circulation cells near the outer boundary in the Navier-Stokes 



fiuid and six in the normal HVBK component. The magnitude of the torque increases 
with the number of circulation cells, because more circulation cells imply steeper radial 
velocity gradients. We observe this in the quantity dNi 2/d{cos6), which measures the 
differential contribution to the torque as a function of colatitude and is plotted in Figure 



17 for the inner and outer spheres. For example, we find |dA^i/(i(cos6')| < 0.1 for the 
Navier-Stokes fiuid, but |diVi/d(cos6')| is as large as 0.5 ( at 9 ^ 75°, 110°) for the 
HVBK superfiuid with HV mutual friction. From equation (3.17), we see that A^i and 
A'2 in the HVBK superfiuid are greater due to larger contributions from the first term 
in (3.171, viz. dvn^/dr. For example, a,t r — R2 and 6 w 45°, we find dv^/dr « —0.73 
for the Navier-Stokes fiuid and {dvn,p/dr) « —3.46 for the normal component of the 
HVBK superfiuid, whereas the second term in (3.171 is similar in both systems, viz. 
Vn^/r « v^/r « 0.49. 

Behaviour of the sort just described was predicted by Henderson & Barenghi (19951, 
who solved the HVBK equations numerically inside infinitely long, differentially rotating 
cylinders. They too observed that, as Re increases, the tension force diminishes, and 
the friction force dominates. Inside the circulation cells of the normal fiuid, the ratio 
|T|/|Fhv| decreases with increasing Re. Henderson & Barenghi (19951 suggested that, 
at higher Re , the mutual friction locks together v„ and Vs , so that the streamlines of the 
superfiuid and normal fiuid are similar. 
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In order to quantify how the two- fluid coupHng changes with Re, consider a HVBK 
superfluid with the same parameters as in Figure [T6| 3 and [16]:, but with i?e = 178 instead 
of Re ~ 10*. Figure 18 compares a sequence of snapshots of the meridional streamlines 
for Re — 178 and lO^t times 10 ^ 60. For Re = 178, the normal component differs 
markedly from the superfluid component, except during the early stages of the evolution 
{t ^ 10). Eventually, at i ^ 30, the normal fluid settles down into a permanent 0- vortex 
state, while the superfluid develops 3-4 vortices near the equator. For Re = 10*, on the 
other hand, the normal and superfluid components display similar flow patterns, with 
the same number of vortices in approximately the same locations. This occurs because 
the HV mutual friction progressively dominates the vortex tension as Re increases and 
also as time passes. Quantitative evidence is presented in Figure 
is plotted as a function of colatitude in the equatorial region 80' 
boundary of the outer sphere (r = i?2)- However, this initial transient soon disappears, 
and the inequality reverses. At t = 10, we find (|T|/|FHv|)i^e=i78 ^ (l'^l/I^Hv|)i^e=io'' 
except at 6* « 85° and 6 w 95° (see Figure [T9|i) . At t = 20 

(|t|/|Fhv 
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where |T|/|Fhv| 
9 < 100°, at the 



< 



we find 
} < 93°. 



|FHv|)i?e. 



:178 



> 



For t ^ 30, we find 



)i?e=io«'' except in the narrow region 87° 
(|T|/|FHv|)i?e=i78 > (|T|/|FHv|)i?e=io* throughout the equatorial region 80° < 6* ^ 
100° (except for a brief reversal at i « 50). Thus, at low Reynolds numbers {Re < 10'^) 
and at early times, the tension force dominates the mutual friction throughout most 
of the fluid, while the opposite is true at high Reynolds numbers and late times. The 
stiffness of the superfluid vortex array, encoded in the tension force, prevents Vs from 
following v„, whereas, when the mutual friction dominates, copies v„ more closely. 

As the tension is less important at higher Re, one expects the superfluid component 
to influence the overall dynamics less as Re increases. This is reflected in the torque. 
For a viscous fluid at Re = 178, the torque is half that for a superfluid at the same 
Reynolds number. For a superfluid with Re — 10*, the torque doubles when compared 
with a viscous fluid with the same Reynolds number. In Section [5. 5| we quantify how the 
boundary condition on the superfluid component (indirectly) affects the torque. 

5.5. Effect of the boundary conditions 

The streamlines of the normal and superfluid components resemble each other ever more 
closely as Re increases, suggesting that the frictional coupling is responsible, as argued in 
Section [5. 4[ Nevertheless, it is important to check how much of the similarity arises from 
imposing the no-slip boundary condition on the superfluid component at r = R2, 
which in turn imposes zero counterflow (v„ = Vg) at r = Ri, i?2. This matters, because 
it can be argued that the no-slip condition is artiflcial. The physically correct boundary 
conditions on are still uncertain, lying somewhere between the following two extremes: 
(i) quantized vortex lines slide freely along the surface, thereby terminating perpendicular 



to it, as expressed by (3.11 1; or (ii) quantized vortices are pinned to the surface, so that 
\g does not slip relative to the surface, as expressed by (3.101 ( Khalatnikov||196"P) . 

To clarify these matters, we repeat two of the no-slip simulations in Section 5.4 (with 
Re — 178 and Re — 10*, as well as J = 0.5 and Afi = 0.3) such that the superfluid satisfies 
the perfect-slip boundary condition (3.101 instead of the no-slip condition (3.11 1. Perfect 



slip implies that the vortex lines terminate perpendicular to the surfaces of the inner and 
outer spheres. The normal fluid satisfies the no- slip condition (3.6 1. Both v„ and are 



initialized to the Stokes solution (3.4 1 



Figure |20] compares the results for no slip and perfect slip at t = 18. For both low 
{Re = 178) and high {Re — 10*) Reynolds numbers, the large-scale structure of the fiow 
is the same under both kinds of boundary conditions. For example, in both Figures |20^ 
(perfect slip) andl20b (no shp), exhibits three large circulation cells: one near the pole, 
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centered near the inner sphere, at « 30 and r « 0.6; one at mid-latitudes, centered 
near the outer sphere, at 6 45° and r k, 0.8; and one at the equator, whose diamater is 



half that of the polar cell. Similar structures are also observed at Re — 10 in Figures 20 



(perfect slip) and 20 i (no slip). On the other hand, the detailed internal structure of the 
cells does depend on the type of boundary condition employed, expecially at lower Re. 
For example, when there is no slip, the centers of the circulation cells develop additional 
small vortices, and the streamlines at r = i?2 become jagged. 

In conclusion, therefore, the global resemblance of the v„ and streamlines at high 
Re found in Section |5.4| is not an artifact of imposing no slip on Vs at the boundaries. 
It is observed equally when perfect slip is allowed. Indeed, the choice of boundary condi- 
tions affects the Vg flow pattern only as far as the small-scale structure in the cell cores 
is concerned. (Of course, the v„ streamlines are almost independent of the boundary 
condition used for the superfluid.) If HV mutual friction is replaced by (weaker) GM 
mutual friction, the resemblance lessens, with v„ tending to look like a Navier-Stokes 
flow (Figure 14 \) and tending to look like an uncoupled superfluid (Figure [14)3 ) . 

The torque on the outer and inner spheres is plotted versus time in Figures [21^ - |21| 3, 
for perfect slip (solid curve) and no slip (dashed curve). A viscous Navier-Stokes fluid 
is also plotted for comparison (dotted curve). The perfect-slip boundary condition on 
roughly halves the amplitude of the torque compared to the no-slip boundary condition. 



6. Nonaxisymmetric spherical Couette flow 

We take advantage of the three-dimensional capabilities of our numerical solver to 
investigate two systems that exhibit nonaxisymmetric flow (requiring spectral resolution 
> 12): (i) a spherical, differentially rotating shell in which the rotation axes of 
the inner and outer spheres are mutually inclined; and (ii) a spherical, differentially 
rotating shell in which the outer sphere precesses freely, while the inner sphere rotates 
uniformly or is at rest. These systems have never been studied before. We use standard 
vortex identification methods, introduced by phong et al. (1990) in viscous flows, to fully 



characterize the three-dimensional vortex structures we encounter — the first time this 
has been done for an HVBK superfluid. One incidental outcome of the work is to conflrm 
that our numerical method can resolve flne structures in superfluid flow. 

6.1. Characterizing the flow topology 

To understand the topology of a three-dimensional flow, one must classify the vortices 
it contains. In its simplest guise, a vortex coincides with a local pressure minimum, 
where the centrifugal acceleration balances the radial pressure gradient. However, this 
criterion fails in situations where unsteady strains or viscous effects (at low Reynolds 
numbers) balance the centrifugal force (Chong et al. 1990 Jeong & Hussain 19951. A 
second simplistic approach is to use |a;| = |V x v| to identify local concentrations of 



vorticity (Metcalfe et al. 19851. However, this method can lead to confusion in wall- 



bounded flows, like ours, where the vorticity produced by a wall-driven background shear 
can mask the main vortical structures of the flow, and one needs to know the location of 
the vortex core beforehand ( Robinsoii||1991 1 . 

Several sophisticated alternatives to the simple pressure minimum and vorticity mag- 
nitude tests have been developed. Most are based on invariants of the velocity gradient 
tensor ( Hunt et a?.||1988" l. We describe two tests below: the discriminant deflnition of a 
vortex, which is insensitive to numerical error, especially in Couette geometries (Frana 



et al. 2005); and the A2 definition, which is suited better to open geometries, e.g. the 
flow past a rotating sphere ( Giacobelio]|2005 ) . 
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6.1.1. Discriminant definition of a vortex 

The velocity gradient tensor Aij = dvi/dxj measured by a nonrotating observer that 
moves locally with the fluid can be decomposed into a symmetric {Sij, or rate-of-strain 
tensor) and an antisymmetric (Wij, or rate-of-rotation tensor) part: 

The eigenvalues A are roots of the characteristic polynomial A"^ + P/iA^ + Q^A + Ra = 0, 
whose coefficients are defined by 

Pa = -tr(A,,) = ~S,,, (6.2) 
Qa = 1[Pa~ tr(4.)] - ^(^1 - - W.,jW,,), (6.3) 

Ra = -det(A,j) = ^(-Pj + SPaQa - 5'jfe^fc» - 3W,jW,kWk^). (6.4) 

These quantities are invariant under any nonrotating coordinate transformation. The 
quantity Pa is the trace of the velocity gradient tensor (that is, the continuity equation); 
in an incompressible fluid, one has Pa = 0. The quantity Qa measures the excess of 
rotation over strain. The sign of Ra governs the stability of the flow (see below). The 
scalar invariants Ra and Qa can be combined to form the discriminant Da = Q\ + 
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• If Da > at some point in the flow, Aij at that point has one real and two complex 



conjugate eigenvalues. Following the nomenclature of |Chong et a?.| ( |1990[ ), we say that 



such points are focal in nature. Streamlines wrap around the axis of the real eigenvector 
and describe a spiral in the plane spanned by the two complex eigenvectors. The sense of 
the spiral is determined by the sign of Ra- For Ra > 0, the trajectories are attracted to- 
wards the axis of the real eigenvector; the point is an unstable focus/contracting (UF/C) 
( Chong et a?. [1990 ^. For Ra < 0, the trajectories are repelled away from the eigenvector; 



the point is a stable focus/stretching (SF/S). 

• If Da < 0, all the eigenvalues are real. We say that such points are strain dom- 



inated (Soria & Cantwell 19941. Streamlines either approach or flee the point along 
the three independent, intersecting, real eigenvectors. When projected onto the three 
planes spanned by the eigenvectors, the trajectories asymptotically approach the eigen- 
vector axes along "parabola-like" or "hyperbola-like" paths, depending on the sign of 
Ra- If Ra < 0, we get a stable node/saddle/saddle (SN/S/S); if Ra > 0, we get an 
unstable/node/saddle/saddle (UN/S/S). The degenerate case Ra — corresponds to a 
two-dimensional flow. 

In this paper, we plot topological isosurfaces according to the following colour scheme. 
Regions with Da > which are SF/S and UF/C are coloured yellow and light blue 
respectively. Regions with Da < which are SN/S/S and UN/S/S are coloured orange 
and green respectively. 

6.1.2. A2 definition of a vortex 



Jeong & Hussain ( 1995 1 realized that the existence of a pressure minimum is not 
a sufficient and necessary condition for a vortex core to be present. For example, an 
unsteady flow can create a pressure minimum even where there is no vortex. On the 
other hand, centrifugal forces can cancel viscous forces perfectly (e.g. Karman's viscous 
flow, or a Stokes flow at low Re), thereby eliminating a pressure minimum even though 



a vortex is present ( Jeong fc Hussain||1995 ) 



The A2 criterion seeks to locate the pressure minimum in a plane perpendicular to the 
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vortex axis by eliminating unsteady strains and viscous stresses from the Navier-Stokes 
equations. Decomposing A^j into symmetric and anti-symmetric parts and substituting 
into the Navier-Stokes equation, we obtain, for the symmetric part. 



dt 



dxk 



1 d^S,, 
Re dxhdxh 



Sk,S,k + Wk,W,k = 



(6.5) 



The Hessian of the pressure, d^p/dxidxj, contains information about pressure minima. 
The existence of a pressure minimum requires the Hessian to have two positive eigenvalues 
( Jeong fc Hussain|1995 1. Ignoring the unsteady strain [first term in (6.5 1] and the viscous 
stress [third term in (6.5)], a vortex is defined as a connected region with two negative 
eigenvalues of the tensor + W^, which is symmetric and has real eigenvalues only. 
Hence, if we call the eigenvalues Ai, A2, and A3 and order them such that Ai ^ A2 ^ A3, 
we define a vortex to be a region where one has A2 < dJeong fc Hussain||1995[ ). 



6.2. Misaligned spheres 

In this section, we consider a spherical shell filled with superfluid {6 — 0.3, Re — 10"^ or 
lO"*^), whose inner and outer surfaces rotate at the same angular speed, Qi — Q2 = 1-0, 
but about different rotation axes. The outer sphere rotates about an axis inclined with 
respect to the z axis, in the x-z plane, by an angle = 3°, while the inner sphere 
rotates about the z axis. We present a gallery of some of the flows excited in this system 



and characterize their topologies employing the methods in Section 6.1 Our results are 
merely a first attempt at this problem; a lot remains to be learnt about the nonlinear 
physics behind such flows. 

No-slip boundary conditions (3.6 1 are imposed on v„, while the superfluid satisfies 



either perfect slip (3.11 1 or no slip (3.101. The mutual friction takes the GM form (2.8 1 



We fail to obtain stable evolution for HV mutual friction (2.6) or large inclination angles 
(^0 ^3°); the simulation becomes unresolved in (j) unless ^ 20 is used, straining our 
computing budget jf] 

6.2.1. Topology of the flow 

Investigations of the flow past a sphere ( Giacobello|2005 1 found that the A2 criterion is 



better at identifying nonaxisymmetric vortex structures than the discriminant criterion. 



validating the findings of Jeong & Hussain ( |1995) , who noticed that vortical structures 
identified using Da tend to have noisy boundaries. Recently, however, investigations of 



Taylor-Gortler vortices in cylindrical containers by Frana et al. ( 2005 1 revealed that the 



vortex contours identified by the A2 criterion can be severely perturbed by numerical 
noise, unlike the discriminant criterion. We find this too. The discriminant criterion is 
better suited to enclosed geometries than the A2 criterion, while the opposite is true for 
open geometries. 

In Figures 22i-c, we plot isosurfaces of A2 (A2 = —0.5, —1, —3) for the normal com- 
ponent in superfluid SCF with Re = 10^, at t = 50. It is hard to discern the regions 
of the flow which are vorticity dominated. Moreover, when varying the isosurface from 
A2 = —0.5 to A2 = —3, we find that its overall shape changes greatly, which is undesir- 
able; the visualization method should be insensitive to the threshold ( Giacobello||2005 1 . 
By contrast, when we use the discriminant criterion, as in Figures |22| i-f, focal regions 
in v„ are cleary visible, and their shape is preserved despite varying the threshold over 
four decades (lO""^ -Da < !)• 



f We note here that post-processing and visualization of the data is a very time-consuming 
task when dealing with nonaxisymmetric flows. 
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The discriminant criterion allows us to diagnose the topology of the inviscid superfluid 
component as well. In Figure 23 we present isosurfaces of Dj^ = 10^'* (Figures 23i-d) 

103 



and Da 



-10- 



(Figures 23 ^-h) for in superfluid SCF with Re — 10'^ and no- 
Throughout most of the volume, the flow is focal, or 



slip boundary conditions on 
vorticity-dominated. Strain-dominated regions, shown in orange, also exist, but are less 
widespread. They have a threaded structure (Figures |23| 3-h) , which is maintained when 
we increase the Reynolds number to Re — lO**, because is weakly coupled to v„ via GM 
mutual friction (see Section 5.4 1. Perfect-slip boundary conditions on Vs do not alter the 
topology. The normal fluid dynamics, on the other hand, is almost completely dominated 
by vorticity, as Figures [24^-d show. Strain-dominated regions are only detected in small 
regions close to the poles (see Figures [24^-h) . This is natural: the differential rotation is 
small (^iCOsGq — fii « 0.01), so the strain applied to the fluid is small. 

The changes in the flow from one snapshot to the next are hard to discern in Figures 
23 and 24 In order to reveal them more clearly, we plot isosurfaces of {yis\ 
at times 10 ^ t ^ 50 in Figure 25 Positive values [(a;s.n)<; 
blue; negative values [(a; 



and (a;„)^ 
0.1] are drawn in light 
)0 = —0.1] are drawn in orange. Figures 25i-e show how the 
normal fluid component evolves. The wedge-shaped isosurface (wn)^ = —0.1 develops 
a pointy extension at the equator that spreads clockwise. Likewise, for the superfluid 
component, the isosurface (ws)^ — —0.1 spreads clockwise in an equatorial band located 
at 60° < < 150°. Note that, although the changes in the flow are easier to see, the 
isosurfaces are threshold-sensitive. For example, the wavy contours in Figure 25i-e are 
not visible for (ws.n)^ = ±10^^. 
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(b) pure one-component inviscid superfluid, (c) normal fluid component of an HVBK superfluid 
with HV friction, (d) superfluid component of an HVBK superfluid with HV friction, (c) normal 
fluid component with GM friction, and (f) superfluid component with GM friction. All snapshots 
are taken at t = 6. SCF parameters: 6 = 0.5 and AQ = 0.3. 
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Figure 15. Torque on the inner (upper plot) and outer (lower plot) spheres, for a viscous fluid 
(solid curve) , a superfluid with HV mutual friction force (dashed curve) , and a superfluid with 

GM mutual friction (dotted curve). In all cases, we take Re = 10*, 5 = 0.5, and AO = 0.3. The 
solid and dotted curves differ by one part in 10'', and are almost indistinguishable on the scale 
of the plot. 
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Figure 16. Meridional streamlines in SCF for (a) viscous Navier-Stokes fluid, (b) normal fluid 

component, with HV mutual friction, and (c) superfluid component with HV mutual friction, 
with Re = 10^. All snapshots are taken at i = 20. SCF parameters: 5 = 0.5 and Af2 = 0.3. 
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Figure 17. Differential torque dNz /d(cosO} as a function of colatitude 6 for the inner (upper 
plot) and outer (lower plot) spheres, for a viscous Navicr Stokes fluid (solid curve), a superfluid 
with HV mutual friction (dashed curve), and a superfluid with GM mutual friction (dotted 
curve), at time t = 20. SCF parameters: Re = 10*, d = 0.5, and An = 0.3. 
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Figure 18. Meridional streamlines in superfluid SCF at low and high Reynolds numbers. Left 
half of figure: Re = 178, at times (a) t = 10, (b) t = 20, (c) t = 30, (d) t = 40, (e) t = 50, and 
(f) t = 60. Right half of figure: Re = 10*, at times (g) t = 10, (h) t = 20, (i) t = 30, (j) t = 40, 
(k) t = 50, and (1) t = 60. In each snapshot, the left-hand and right-hand quadrants display 
the normal and superfiuid components respectively. SCF parameters: S = 0.5, Afi = 0.3, and 
Re, = 10^ 
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Figure 19. Ratio of the tension force (T) to the HV mutual friction force (Fhv) as a function 
of colatitude 9 (in °), at the boundary of the outer sphere, at times (a) t — 10, (b) t = 20, (c) 
t = 30, (d) t = 40, (e) t = 50, and (f) t = 60. SCF parameters: 5 = 0.5, AO = 0.3, Re = 178 
(solid curve), and Re = 10^ (dashed curve). Generally speaking, |T|/|Fhv| decreases with Re 
and t. 
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Figure 20. Meridional streamlines in superfluid SCF for the normal (left panels) and superfluid 
(right panels) components, at time t — 18, for S = 0.5, AQ — 0.3, and Rcs = 10''. The panels 
illustrate the effect of boundary conditions, (a) Re = 178. perfect-slip boundary condition (3.101 
on Vs ; (b) Re = 178, no-slip boundary condition ( 3. 1 1 1 on Vs ; (c) Re = 10'* , perfect-slip boundary 
condition (3.101 ; and (d) Re = 10'*, no-slip boundary condition. 
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time time 

Figure 21. Evolution of the inner torque (upper plots) and outer torque (lower plots) for super- 
fluid SCF with <5 = 0.5, AQ = 0.3, and Reynolds number (a) Re = 178 and (b) Re = lO". The 
solid, dashed, and dotted curves correspond to an HVBK superfluid with perfect-slip boundary 
condition on Vs, an HVBK superfluid with no-slip boundary condition on Vs, and a viscous 
Navier-Stokes fluid, respectively. 
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(a) 





(b) 





(c) 



(f) 





Figure 22. Instantaneous flow topology for the normal component in superfluid SCF, with 
Re = 10'^, S = 0.3, Af2 = 0, and = 3°. Isosurfaces in orange for (a) A2 — —0.5, (b) A2 = —1, 
and A2 = "3, and in light blue for (d) Da = 10"*, (e) Da = 10"\ and (f) Da = 1. The Da 
criterion is to be preferred over the A2 criterion as it is less sensitive to the threshold. Spectral 
resolution and fllter parameters: {Nr,Ng,N^) = (100,200,12), (7r,7e) = (10,10). 
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Figure 23. Instantaneous flow topology of the superfluid component in superfluid SCF with 
Re = 10^, 5 — 0.3, An = 0, So = 3°, and no-shp boundary conditions on Vs. Isosurfaces in Hght 
blue for Da = 10"* at (a) t = 20, (b) t = 30, (c) t = 40, and (d) t = 50; and in orange for 
Da = -10"^ at (e) t = 20, (f) t = 30, (g) t = 40, and (h) t = 50. Spectral resolution and fllter 
parameters: (iV,., iVe, iV^) = (100,200,12), (7r,7e) = (10,10). 
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Figure 24. Instantaneous flow topology of the normal fluid component in superfluid SCF with 
Re = 10*, 5 — 0.3, AQ, = 0, do = 3°, and no-slip boundary conditions on Vs. Isosurfaces in light 
blue for Da = 10~* at (a) t = 20, (b) t = 30, (c) t = 40, and (d) t = 50; and in orange for 
Da = -10"* at (e) t = 20, (f) t = 30, (g) t = 40, and (h) t = 50. Spectral resolution and fllter 
parameters: {Nr,Ne,N^) = (100,200,12), (7r,7e) = (10,10). 
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Figure 26. Inner (left) and outer (right) torque versus time in superfluid SCF with no-slip 
(solid curve) and perfect-slip (dashed curve) boundary conditions on V3, with Re = 10^, 5 = 0.3, 
ASl = 0, ^0 = 3°, and GM mutual friction. From top to bottom, the plots show the (a) x, (b) y, 
and (c) z components of the inner torque, and the (d) x, (e) y, and (f) z components of the outer 
torque. Spectral resolution: [Nr,Ng,N^) = (100,200, 12). Filter parameters: (7r,7e) = (10, 10). 



6.2.2. Unsteady torque 



We present the evolution of the torque on the outer sphere in Figures 26 {Re — 10 ) 
and[27](i?e = 10^). For Re = 10^, we have - lOTV^, and - ^QN^, as expected for 
9o. Moreover, tends to a constant value « 1.7 x 10"'^ for the outer torque 



sma 

and Nz ~ lAx 10~^ for the inner torque at t ^ 20. The boundary condition on Vg has a 
negligible effect. When it is changed from no slip to perfect slip, Nz decreases by < 0.3 
%, and Ny decreases by 1 % (dashed curve in Figure 26 d). For Re = 10^, the torque 
tends to a constant value more gradually than for Re = 10"^. The differences between 
no slip (dashed curves) and perfect slip (solid curves) are slightly greater; and Ny 
(see Figures 27 i-b and [27}i-e) are ^ 2 % larger for no slip. Again, the dominant torque 
component is {> Ny > N^). 

The differences in the torque components arise from asymmetries in the flow. In an 
axisymmetric flow, the greatest contribution to the torque comes from regions containing 
a larger number of tightly packed circulation cells. The same is true in a nonaxisymmetric 
flow. We calculate the torque from N = J dS{Trgee + Tr^eif,), where dS denotes the area 
element on the sphere and r^j are the shear stresses. We have Tro 7^ and r^^ 7^ in a 
nonaxisymmetric flow, giving N^ ^ and Ny ^ 0. Note that N^ ^ Ny when avergared 
over time, since the rotation axis is in the x-z plane. 




Figure 27. Inner (left) and outer (right) torque versus time in superfluid SCF with no-shp 
(sohd curve) and perfect-shp (dashed curve) boundary conditions on Vg, with Re = 10*, 5 = 0.3, 
ASl = 0, ^0 = 3°, and GM mutual friction. From top to bottom, the plots show the (a) x, (b) y, 
and (c) z components of the inner torque, and the (d) x, (e) y, and (f) z components of the outer 
torque. Spectral resolution: (iV^, iVe, iV</,) = (100,200, 12). Filter parameters: (7r,7e) = (10, 10). 



6.3. Free precession 

In this section, we consider a spherical rotating shell filled with superfluid, where the 
outer sphere precesses freely, while the inner sphere rotates uniformly. We exaggerate 
the biaxiality of the outer shell, taking Hp ~ 1.0 for the body-frame precessional fre- 
quency (defined in Section 3.5.41 and fi' = 2.0 for the inertial-frame precession frequency 
(Landau & Lifshitz|[l969 Jones fc Andersson|2001 1. This allows us to investigate all the 



time-scales comprising the precession dynamics using simulations of reasonable duration, 
something that would be impossible for Up <C O'. The fixed angular momentum vector of 
the outer sphere points in the z direction in the inertial frame of the inner sphere (which 
rotates with fii — 1.0). A n expression for the velocity of every point on the outer sphere 
is given in Section 3.5.4 We consider a relatively low Reynolds number. Re = 10'^, with 
6 = 0.3 and no-slip boundary conditions on v^. 

6.3.1. Topology of the flow 

The topology of the flow is illustrated in Figure [28] for the normal fluid component and 
in Figure [29] for the superfluid component. Unlike the misaligned spheres in Section [6?2] 
this flow is influenced equally by strain and vorticity. The UF/C topology is slightly more 



prevalent (see the light blue isosurfaces in Figures 28i-d) than the SN/S/S topology 
in the normal fluid component. In the superfluid component, the UF/C and SN/S/S 
topologies are equally prevalent. The UF/C regions (in light blue in Figures 28i-d) 
exhibit a comphcated brick-like structure (for = 12), while the SN/S/S regions are 
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Figure 28. Instantaneous flow topology for a rotating superfluid contained within a freely 
precessing outer sphere and uniformly rotating inner sphere. Discriminant isosurfaces are shown 
for the normal fluid component, for Da = 10~* (light blue) at (a) t — 20, (b) t — 30, (c) t = 40, 
and (d) t = 50, and for Da = -10"* (orange), at (e) t = 20, (f) t = 30, (g) t = 40, and (h) 
t = 50. The mutual friction is of GM form. Simulation parameters: Re = 10'^, S — 0.3, fii = 1.0, 
Q' = 2.0, Qp — 2.0, = 3°, and no-slip boundary conditions on v^. Spectral resolution and 
fiher parameters: {Nr,Ng,N^) = (100,200,12), (7r,7fl) = (8,8). 



more filamentary. The superfluid component is similar to the normal fluid component but 
has smoother isosurfaces (see Figure 29 1 , so it is doubly difficult to distinguish transients 
in the flow. 

When we plot isosurfaces of vorticity, in the same manner as in Figure 25 the results 
are unsatisfactory. The positive and negative isosurfaces are tightly interleaved and it 
is hard to make out the underlying topology. However, the results improve dramatically 
when we subtract the vorticity of the Stokes solution from the total vorticity and project 
Au)n,s along the instantaneous principal axis of inertia, e3(i), of the outer sphere (defined 
in Section 3.5.4| . We present isosurfaces for Au;„ -e^ = ±0.1 in Figures 30i-d; as before, 
positive (negative) isosurfaces are coloured light blue (orange). Similarly, we present 
isosurfaces for Aa;, 



• 63 = ±0.1 in Figures 30' 



-h. We observe that isosurfaces of AuJn ■ ^3 
form two interlocking ribbons of opposite sign which attach {t = 20), detach {t = 40), 
and attach again (t — 50) at two equatorial points (one of which is framed by a black 
circle). In contrast, isosurfaces of -63 exhibit a tongue- like structure in the equatorial 
plane (framed by a black square) , which grows clockwise from t — 20 to t — AO and finally 
develops sawteeth at i = 50 (see Figures [SOjj-h). We suspect that these three-dimensional 
structures are not completely developed by t = 50, i.e. we are observing transients, but 
computational limitations prevented us from extending the runs at the time of writing. 




Figure 29. Instantaneous flow topology for a rotating superfluid contained within a freely 
precessing outer sphere and uniformly rotating inner sphere. Discriminant isosurfaces are shown 
for the superfluid component, for Da = 10"* (light blue), at (a) t = 20, (b) t = 30, (c) t = 40, 
and (d) t = 50, and for Da = -10"* (orange) at (e) t = 20, (f) t = 30, (g) t = 40, and (h) 
t = 50. The mutual friction is of GM form. Simulation parameters: Re = 10 , S — 0.3, Qi = 1.0 
Q' = 2.0, Qp — 1.0, 6p = 3°, and no-slip boundary conditions on v^. Spectral resolution and 
filter parameters: {Nr,Ng,N^) = (100,200,12), (7r,7e) = (8,8). 
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6.3.2. Unsteady torque 



In Figure 31 we plot the viscous torque exerted by the fluid on the inner (Figures 
31i-c) and outer (Figures 31i-f) spheres for Re = 10^. On the inner sphere, we find 



~ lO^iVj, lO'^Ny. On the outer sphere, we find ^ ^ lONy. The outer torque 
increases linearly up to t « 20, when it reaches \Nz\ « 0.03, while the inner torque 
decreases to \Nz\ « 0.01 over the same interval and oscillates persistently (AjTV^j « 10~^, 
period w 3). The other torque components oscillate persistently for 2. For example, 
(N2)x has constant amplitude (fa 4 x 10~^) and period (« 3), (Ni)^; has a smaller 
amplitude (« 2 x 10"'*) but the same period, (Ni)^ has period w 3 and amplitude 
ranging from 2 x 10^"* to 4 x 10^"*, and (N2)y has amplitude ranging from 2 x 10^"^ to 
5 X 10-*. 



7. Conclusion 

Superfluid SCF, like its classical (Navier-Stokes) counterpart, is controlled by three 
global parameters: the dimensionless gap width 6, the Reynolds number Re, and the rota- 
tional shear Afi/fi. In addition, it is a function of the form (isotropic versus anisotropic) 
and dimensionless amplitude of the mutual friction and vortex tension forces. In this 
paper, we solve numerically the HVBK equations describing superfluid SCF for a range 
of S, Re, and Afi/n and study the time-dependent behaviour of the resulting flow. The 
numerical solver is based upon a pseudospectral collocation method. Special attention is 
paid to the pole parity problem and to controlling the growth of global oscillations (due 
to the Gibbs phenomenon) by flltering out high spatial frequencies spectrally. The solver 
accurately resolves flows covering the parameter range 10 ^ i?e ^ 10^, 0.2 ^ (S ^ 0.9, and 
^ An/n < 0.3. Grids with resolution {Nr,Ne,N^) = (150,400,4) and (100,200,12) 
are adopted for the most challenging problems we attempt in two and three dimensions 
respectively. 

In two dimensional superfluid SCF, persistent quasiperiodic oscillations are always 
observed in the torque during steady differential rotation (after initial transients die 
away), with typical period ~ and fractional amplitude ^ 10~^. The oscillation 

amplitude increases as Re increases. The viscous torques exerted by a Navier-Stokes fluid 
and an HVBK superfluid with GM mutual friction differ by 6 %. However, the torque 
roughly triples for HV mutual friction. The meridional streamlines are more complicated 
for HV friction, with more small circulation cells near the outer sphere, and the amplitude 
of the torque oscillations is greater. 

In three dimensional superfluid SCF, nonaxisymmetric vortex structures are classifled 
according to topological invariants. The discriminant criterion is more instructive than 
the A2 criterion. For misaligned spheres, the flow is focal (vorticity-dominated) through- 
out most of its volume, except for thread-like zones where it is strain-dominated near the 
equator (inviscid component) and poles (viscous component). A wedge-shaped isosur- 
face of vorticity rotates around the equator at roughly the rotation period. For a freely 
precessing outer sphere, the flow is equally strain- and vorticity-dominated throughout 
its volume. Unstable focus/contracting points are slightly more common than stable 
node/saddle/saddle points in the viscous component but not in the inviscid component. 
Isosurfaces of positive and negative vorticity form interlocking poloidal ribbons (viscous 
component) or toroidal tongues (inviscid component) which attach and detach at roughly 
the rotation period. Persistent torque oscillations are observed in all the three dimen- 
sional flows considered, with period ~ 6^1^. 

A detailed knowledge of the global superfluid hydrodynamics inside a neutron star is 
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Figure 31. Evolution of the viscous torque on the inner (left) and outer (right) spheres for a 
rotating superfluid contained within a freely precessing outer sphere and uniformly rotating inner 
sphere. From top to bottom, the plots show the (a) x, (b) y, and (c) z components of the inner 
torque and the (d) x, (e) y, and (f) z components of the outer torque. Simulation parameters: 
Re = 10"^, 5 = 0.3, Q! = 2.0, f2p = 1.0, fii — 1.0, Op — 3°, no-slip boundary conditions on Vs, 
and GM mutual friction. Spectral resolution and filter parameters: {Nr, No, N^) — (100, 200, 12), 
(7r,7e) = (8,8). 



needed to understand the origin of the timing irregularities — glitches and timing noise 



observed in over 100 radio pulsars to date (D'Alessandro et ar||1995 Shemar & Lyne 



1996 



2005 



Lyne e< a?.|[2000{ |Hobbs|[2002{ [Scott et aL||2003| |Hobbs et a?.||2004[ [Peralta et al. 



2006a|& Melatos et al. 20071. Glitches are characterized by a s udden increase in 
the angular velocity of the pulsar, in the range 10~^^ < Afi/fi < 10~^ (Lyne & Graham- 
Smith|2006| |Peralta|2006| . Pulsars also exhibit non-Gaussian fluctuations in pulse arrival 



times over many years, known as timing noise. The long relaxation time following a glitch, 
and the temperatures in a neutron star, independently imply that the interior of the star 
is a neutron superfluid ( Lorimer fc Kramer||2004 Lyne fc Graham-Smith|[2006 ) . 

Recently, it was shown that the global pattern of superfluid circulation in a neutron 
star exerts a dramatic influence on its rotation and may play a central role in explain- 
ing the phenomena of glitches and timing noise (Peralta et al. 2005 2006a|& I. For this 
reason, it is important to understand more fully the axisymmetric and nonaxisymmetric 
dynamics of superfluid SCF. The results in this paper represent a flrst step along this 
path. One particular consequence is that, if the meridional circulation is fast enough, a 
vortex tangle {superfluid turbulence) is alternatively created and destroyed in the outer 
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core of the star (and indeed any spherical container). For example, before a glitch, differ- 
ential rotation in the outer core drives a nonzero, poloidal counterflow which excites the 
Donnelly-Glaberson instability (DGI) ( Glaberson et aL|1974 Swanson et aL|1983|[Tsub^ 
ota et aL|2003 1, and the vortices evolve into an isotropic tangle of reconnecting loops. In 
this regime, the friction force is of GM form, coupling the normal and superfluid com- 
ponents loosely. Right after a glitch, the differential rotation ceases, so does the poloidal 
counterflow, the vortex tangle decays, a rectilinear vortex array forms, and the mutual 
friction changes to HV form, suddenly locking the normal and superfluid components 
together. In previous simulations done for neutron star parameters ( Peralta et aL||2"005 l, 
we make an order-of-magnitude estimate of the ratio |Fhv|/Fgm| as follows. Ignoring 
B', we find |Fhv|/Fgm| ^ [Bujs{vns - Vs/ R2)]/[A' PnPsvlJ np^]. Taking v^s ^ R2^2, 
with R2 - 10^ cm, ^2 ^ 10^ Hz, and A' - 10"^ we get |Fhv|/Fgm| ^ 10^. From the 
numerical simulations, we obtain a typical value |FhvI/Fgm| ^ lO^i which is similar. 
Nevertheless, we note that the microphysics of the GM force in superfluid turbulence has 
not been worked out fully yet (Jou & Mongiovi l[2004l ). 

The results on superfluid turbulence summarised above are also relevant to laboratory 



experiments by |Tsakadze fc Tsakadze| ( |1972[|l973[|1974[|1975[|1986p , the only systematic 
experimental study of spherical Couette flow in superfluid helium undertaken to date. 



Tsakadze & Tsakadze ( 1972 ) studied the deceleration of axisymmetric vessels made of 
glass and plastic and flUed with ''He (in the temperature range 1.4K ^ T ^ 2.0 K), after 
an impulsive acceleration. They observed "jerky" behavior, reminiscent of pulsar glitches, 
and developed an empirical formula for the relaxation time as a function of the initial 
angular velocity, the normal fluid fraction, and the radius of the vessel. The results agree 
qualitatively with glitch data from the Crab and Vela. These experiments lend support 
to the idea that the neutron superfluid inside pulsars plays an important role in the glitch 
phenomenon. More broadly, however, the Tsakadze experiments — and by extension, the 
theoretical results in this paper — are of general interest in understanding the physics 
of superfluid turbulence in rotating systems ( Barenghi et a?.|2001 l. The spin- up problem 
in helium II, and superfluids in general, is far more complicated than in a viscous fluid, 
because the normal fluid component interacts nonlinearly with the quantized Feynman- 
Onsager vortices in the superfluid component. One interesting effect is that sudden, 
"glitch-like" , spin-up events and other rotational irregularities are associated with patchy 
mutual friction: the DGI is excited in parts of the superfluid (e.g. near the walls, on the 
rotation axis, and at the equator) but not elsewhere (Peralta et al. 2006&). This new 
phenomenon will be investigated further in a forthcoming paper ( Peralta et aL|p007 l. 
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Appendix A. Pseudospectral collocation grids 

The radial direction (r) is discretized using Chebyshev polynomials and a collocation 
scheme ( Boyd|[200T |. The Gauss-Lobatto collocation points in r are defined as (Canuto 
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erZ][l988 l 



nd — 1) 

Xi = -cos Y ' with 1 i iVr, (Al) 

where Nr is the number of radial collocation points. The computational points in 
physical space are related to the Chebyshev grid Xi through the mapping 



R2 ^ Ri 



Ri + R2 



(A 2) 



The angular directions 9 and are discretized using a periodic grid over the intervals 
^ ^ TT and Q ^ (j) respectively. In the azimuthal direction, the collocation points 
are defined as 

_ 27r(fc- 1) 



with 1 fc ^ iVw,, 



(A3) 



where is the number of grid points in (f). In the polar direction, the collocation points 
are defined as 

ttQ- - 1/2) 

^ ' - No ' 
where Ng is the number of grid points in 9. 



with 1 i^j ^ Ne 



(A 4) 



Appendix B. Spectral expansions 

A scalar field s (such as pressure) is expanded as|f] 

( N,.-lNe-l 

J2 J2 J2 aijkTi{r)cos{j9)e''"^ for fc even 

/=0 3=0 k = -N^/2 

EE E o^ijkTi{r)sm{j9)e''"^ for fc odd 

;=0 J=l fc=-Afj./2 



S = < 



(Bl) 



The radial component of a vector is continuous across the poles, but the tangential 
and azimuthal components change sign. Hence the radial, tangential, and azimuthal 
components of the velocity field are expanded as 



Ur = < 



f Nr-lNa-1 N^/2 

E E E PijkTi{r)cosU9)e'''^ for fc even 

/=0 3=0 k=-N^/2 

EE E f3ijkTi{r)sin{j9)e"'* for fc odd 

1=0 3 = 1 k=~N^/2 



(B2) 



(B3) 



( N^-^l Ne N^/2 

EE E lijkTiir)smij9)e''"^ for fc even 

_ . 1=0 j=l k=~N^/2 
Ue - \ M^-lNe-l N^/2 

E E E li3kTi{T)cos{j9)e'''^ for fc odd 

/=0 j=0 k=^N^/2 

f Note that, in ( |A 1[ )-( [A^ , the grid indices run from 1 to Nr^e,ij>, whereas, in ( |B 1[ |-( [B^ , the 
wavenumber indices I, j, k run over different ranges. In the solver, the Navier-Stokes equations 
are discretized according to the prescription in Appendix [C| 



U0 = < 



(B4) 
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( Nr~l Na 

J2 5i,kTi{r) s\nme'^^ for fc even 

1=0 j = l k=-N^/2 

J2 5i,kTi{r)cos{j0)e'''^ for fc odd 

1=0 ]=0 k=-N^/2 

where Ti is the ^-th Chebyshev polynomial, j and k are the and (j) wavenumbers, and 
aijk, f3ijk, li]k and 5ijk are real coefficients. 



Appendix C. Solution Algorithm 

C.l. Numerical differentiation 

Spatial differentiation in r and 9 is carried out in physical space. The first-order r deriva- 
tive is calculated by performing the operation 



9/ 
dr 



with l^i^Nr, 



(CI) 



in which f{xj) is a discrete vector array, and 'v'f^ is the differentiation operator for the 
variable r, defined below. A similar formula apphes for df /dO. 



For the Gauss-Lobatto distribution of radial points defined in Section 3.1 the radial 
derivative operator is defined as 



(-1) 



i+i 



2 sin 



2{Nr~l){i + j -2) 



sm 



T^ji - j) 
2{Nr - 1) 



cos[7rj/(iV^ - 1)] 

2W['Kj/{Nr-l)] 
2{Nr -1)^+1 

6 

2iNr -1)2 + 1 



for i j 



for i <— j < Nr 



for i — j — 1 



(C2) 



6 



for i 



Nr 



The top-half entries of the matrix (1 ^ i,j ^ Nr/2) are more accurately represented 



than the lower-half ones, for a given machine precision e ( Don fc Solomonoff|1995 1. A sig- 
nificant improvement in accuracy is obtained by using the property V^j = —V 
[for {Nr/2) + 1 ^ Nr], and calculating only the more accurate top-half part of ( C 2 1 



This reduces the overall round-off error from 0{Nf.e) to 0{N^e) (Don & Solomonoff 
|1995[ ). Radial derivatives of order n are computed by applying ( |C 2[ ) n times. 

For a periodic grid with 2Ng points, the first-order 6 differentiation operator can be 
expressed as ( Canuto et aZ.||198"8 l 



■^ij — 



l(_l)-.cot^fc^ 
2 ^ ' 2Nf) 



for i =/: j 



for i — j, 



(C3) 
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and the second-order differentiation operator as 



1 



(-l)'-J+isi 



T^{i- j) 
2Ng 



for i ^ j 



(C4) 



12 



for i 



J- 



The form of Dfj depends on the parity of the function f{xj) on which it acts, as explained 
in Section |b] First-order cosine (C) and sine (S) operators can be constructed from (C3l 
as 



i,2Ne + l-] 



for i^Ng 
for l^i.j^Ne, 
and second-order cosine and sine operators can be defined as 



'^ij ~ -^ij ■^i,2Ng + l-j 



i,2Ne+l-] 



'^ij ~ •' ij •^i,2Ng + l-j 



for 1 ^ Ng 

for l^i.j^Ng. 



(C5) 
(C6) 

(C7) 
(C8) 



The operators (C5l-(C8l are apphed in spectral space using a fast Fourier transform 



(FFT). Since the expansions (B 1 |-(B4| are different for even or odd k, the 9 differenti- 



ation is performed separately for even and odd k. For even fc, a C operator is applied; for 
odd k, an S operator is applied . The differentiation operator reverses the parity when 
applied an odd number of times and leaves the parity unchanged when applied an even 
number of times, which needs to be taken into account when computing higher order 
derivatives. 

For the (f> derivatives, the periodicity of the grid makes it easier to perform the dif- 
ferentiation in spectral space. Differentiation in reduces to multiplying each Fourier 
coefficient FFTcj,[f{xk)] = f{xk) by i = V^l times the corresponding wavenumber /c. 



for 1 s$ j < Nr. 



(C9) 



The function f{xj) is transformed back to physical space using an inverse FFT. In order 



to compute (j) derivatives of order n, the operation (C9l is performed n times 



C.2. Temporal discretization 

The most common and efficient method to solve the Navier-Stokes equation in terms 
of the primitive variables (i.e., the velocity and pressure fields) is the fractional step 
approach (Chorin 1968 Brown et al. 20011. This method proceeds in two time steps, 



decoupling the calculation of the velocity and pressure fields. In the first step, an inter- 
mediate velocity field v* ^ is computed from the velocity field at the 7i-th time level, vJJ ^, 
using the momentum equations (2.1 1 and (2.2 1 and ignoring the pressure and the incom- 



pressibility constraint (2.3 1. In the second step, the pressure is calculated by solving a 



Poisson equation involving the intermediate velocity field to obtain the final, divergence- 
free velocity field at the n + 1-th time level, vJJ;+^. The second step can be seen as a 
projection of v„ ,j onto a space of divergence-free vectors. 

Stable evolution is achieved by treating all the viscous terms implicitly, to avoid the 
stringent upper bound on At imposed by the CFL condition ( Voigt et aL||1984" ). Nonlin- 



ear terms, on the other hand, are treated explicitly (Orszag 



1971a I. Upon discretizing 



equations (2.1) and (2.2 1 in time using an explicit third-order Adams-Bashforth (ABB) 



method for the nonlinear terms and an implicit second-order Crank-Nicolson (CN2) for 
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the viscous terms ( Boyd|pOOT I , and advancing the solution from the time level n to the 
time level we obtain 



At 



1 

12 

PsF 
P 



23 



Ps^UJs PsF 

v„ • Vv„ H Da 

plies P 



- 16 



v„ • Vv„ 



Ps'^^s 
pRCs 



+ 



na 



1 



2Re 

J.2 Qj. 



sin^ I 



dr 



sm 



pRes p 

d'^v'' 1 d 
9 7-2 9r 



e, 

9 S-y* 



ri-2 



9r 



(CIO) 



At 



1 

12 



23 



PnF 
P 



v., • Vv., 



T 



PsVlOs 

pRcs 



Vv 



P 

PsVuo, 



T 



- 16 

P 



v., • Vv., 



T 

Res 



Ps'S/uJs 

pRcs 



(Cll) 



where — (e^, eg, e^) is the triad of spherical polar unit vectors. The term Da is defined 
as 



Dr 

Dg 
Da 



d 





sin ^ 






1 


d 




sin^ 


W9 




1 


d 




sin t 


W9 



smi 



sm ( 



sm( 



dVnr^ 




89 ) 




dvne \ 


2 8Vnr 


89 ) 


r2 89 






89 I 


r2 sin'^ ( 



2 dVne 2Vn9COt9 



dVn4, 



r2 89 



2 COS 6* 8vn,j, 



r2 sin 9 



dun 



r'^ sin^ I 



2COS0 dVnB 



(C12) 
(C13) 
(C14) 



Note that the Laplacian operator V2v splits into two parts: 9 derivatives are computed 
explicitly with an AB3 scheme, while r and (j) derivatives are computed implicitly with a 
CN2 scheme. This procedure avoids numerical instabilities ( Bagchi|[2002 1 and allows us 
to rewrite equation (C 10 1 at the time level 7k- as a Helmholtz equation for v* (see Section 

as). 

The previous advection-difFusion step is followed by a pressure correction step based 



on 



,n+l 



At 



n+l 



(C15) 



By taking the divergence of this equation, and imposing the continuity constraint V 
v^"^'^ = 0, we obtain a Poisson equation for the pressure: 



At 



2^n+l 



(C16) 



This equation is solved implicitly for p''^\ at every time step, in order to get divergence- 
free normal and superfluid velocity fields accurate to second order in time. 



A boundary condition is needed to solve equation (C16). Since physical boundary 



conditions do not apply to v* ^ and p^\^, a number of empirical choices have been 



proposed (Kim & Moin 1985 Bell et al. 19891. Unfortunately, all the choices create 



a spurious numerical boundary layer at the edge of the computational domain. The 
simplest and most common choice for the pressure is a homogeneous boundary condi- 
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tion {dplt}sTldf)\T=R^,R.2 

ity field satisfies no penetration, viz. • v, 
for the tangential components can be obtained by rearranging equation ( |C 15| to read 



(Orszag et al. 19861. The radial component of the veloc- 
= • v"+-'^ = 0. A boundary condition 



AWp: 



n+l) 
n,s ) 



lar 



However, Vp"+^ is unknown at this stage of the 



calculation: we wish to solve for it in order to pressure-correct the fields. Streett fc Hus 



sami 



( 1991 1 argue that, since the time-stepping scheme is second-order accurate in time, 



a second-order accurate (or better) estimate for the velocity fields is sufficient to pre- 
serve global second-order accuracy. Expanding VpJJ^^ in a Taylor series about t = 
and approximating dp^^/dt with a first-order backward difference formula, we obtain 
Vp;^+i = iVp^ ,. - Vpn~s^ + 0{At^) and the boundary condition for the 9 and </) compo- 
nents of the ★ velocity fields then becomes 



- At(2Vp" - Vp"+i)] + ©(At^), 



(C17) 



with X = eg or e^. This gives a small slip velocity at the boundary, at the time level n+l, 
of order O(At^). The particular choices for vJJ^^ are described in more detail below. 

Once p""''^ has been calculated, the updated divergence-free velocity fields v"+^ are 
calculated from (C 15). Note that, since there is no diffusive term in equation (2.2), the 



explicit AB3 method is used to advance Vs from the time level n to the time level *. 

C.3. Advancing the solution in time 
Starting from an initial choice of v„ and that satisfies the continuity equations 



(2.3), the time-stepping procedure consists of an advection step (or, more precisely, an 
advection-diffusion step for the normal fluid), in which the linear terms are handled im- 
plicitly (Crank-Nicolson) and the nonlinear terms explicitly ( Adams-Bashforth) , followed 
by a pressure-correction step. Given the solution at the n-th time level, and appropriate 
boundary conditions, the algorithm proceeds as follows. 

(a) We calculate the Fourier and Chebyshev grids, the differentiation matrices, the 
matrix operations that only depend on the spatial grid needed by the Helmholtz and 
Poisson solvers (see below), any exponential filters to be applied to v„ and v^, and the 
pole filter and anti-aliasing filter matrices (see Section 3.4). 



(&) We calculate the nonlinear and linear parts of (C 10 1 and (C 11 1. The time-stepping 



loop starts. The AB3 method initializes v,\ ^ and v,^ j, using a lower-order, Euler method. 

(c) We perform the advection-diffusion step. All the velocity and pressure variables 
are Fourier expanded in (j). A Helmholtz equation (equation C 18 below) is solved in order 



to get the intermediate v* at the time level The AB3 algorithm is used to step forward 
Vs explicitly to the time level -k. The boundary conditions are applied. 

(rf) We a pply filters to the expansion coefficients and a convective filter, as described in 
Section 3.4 We then pressure-correct v* and v* in order to get vJJ+^ and v""*"^, satisfying 



the continuity equation (2.3). Poisson's equation (C 15) is solved for v* and v* and the 



solution is advanced to the time level n -I- 1 using (C 15) 



(e) We take the inverse Fourier transform of the velocity and pressure fields and write 
out their values on the coordinate grid as a restart file (if desired) . 

We now describe in detail steps (c) and (d). The r, and (/) components of the vector 
Laplacian V^v are coupled and a completely implicit treatment of the diffusive terms is 
cumbersome. To get around this issue, in step (d), the linear diffusive terms in equation 



(2.1) are treated semi-implicitly: only the linear terms in the Laplacian are included in 



the implicit temporal discretization (using the CN algorithm) (Bagchi 2002 1 , while the 
nonlinear terms in (2.1) are treated in an explicit way (using an AB3 algorithm). Since 



rAQ > Ar, the 6 components of the Laplacian can be treated explicitly, without affecting 



the stability of the time-stepping algorithm (see Section C.2 1. Moving the nonlinear terms 
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and 9 derivatives of the Laplacian to the right-hand side of equation ( C 10 1 , and the time 



derivatives to the left-hand side, and Fourier transforming all the variables, we arrive at 



a Helmholtz equation for the normal velocity field v„ 
★ which, for the radial component, takes the form 



ivnr,Vne,Vn4>) at the time step 



sm 



dr 



sin^ 9Re 



kd 



dr I \ At '2 
where all Fourier-transformed quantities are indicated by a 



(C18) 



Similar equations can be 



written for v^^g and v^i^- The right-hand side (RHS)^ contains all the non-linear terms 
and 9 derivatives of the Laplacian evaluated at the time step n. This equation is solved 
in spectral space for each wave number fc^,, as described below. 



The absence of a viscous term in equation (2.2 1 allows us to use an explicit method 
for the evolution in time, as indicated in equation (Clll. For = {vsr,Vse,Vs,p), the 



solution advances in time from rt to 7k- according to an AB3 scheme. 



At 



(NLS)^ 



A . n-1 



12 



(C19) 



where (NLS)^ contains all the nonlinear terms in (2.2) 



The velocity fields obtained from equations (C18) and (C19l are not divergence-free, 



and therefore do not satisty the continuity equation. A correction must be made by 



solving the Poisson equation (C 16l for the pressure, which can be written as 




89 



sm 



n.s 

27r 




(C20) 



The corrected velocity fields vJJ+^ and v"+^ are then calculated using (C 15 1, which gives 

l.-WKt'^t- (C21) 



-n+l 



Equations (C 18 1 and (C 20 1 are solved using a matrix diagonalization method (Canuto 



et al. 1988 Trefethen 20011. The Helmholtz and Poisson equations can be written in 



matrix form, viz. 



[A] 



= [RHS\, 



(C22) 



where [w] is a A^^ x array containing the velocity or pressure fields, and [A\ is a 
Nr X Nr m atrix that represents the discrete operators on the left-hand sides of (C18l 
and (C20l. The matrix [B]^ is the transpose of [B], a, Ng x Ng matrix which is zero for 
the Poisson equation (C20l and contains the second and third terms in the Helmholtz 



equation (C 18l. Finally, we have a = 2Re/At for (C 18 1 and a = for (C 20) 



In order to solve equation (C22l by matrix diagonalization ( [Canuto et al. 1988), we 
decompose [A] and [B] in eigenvectors and eigenvalues for the operators Vf. and Vg as 

[A] = [M] [Xr] [M]-' , [B] = [N] [Xg] [Nr' , (C 23) 

where [M] and [A^] are matrices formed from the eigenvectors of [V^] and [Vg] , and [A, ] 
and [Xg] are diagonal matrices formed from the eigenvalues of [V^] and [Vg] , respectively. 



Substituting equation ( C 23 1 in equation ( C 22 1 leads to 

[A, 



M [Xe 



a \u\ 



(C24) 



with [u] = [M]~' [v] [N] and [s] = [M]^' [RHS] [N]. Since [A^] and [Xg] are diagonal 
matrices, this expression can be simplified to [u] — [Xr + Xg — a]~^[s] and written explicitly 
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as 



for 1 i, j < Nrfi. 



(C25) 



Finally, \v] can be obtained from \G2b\ using [u] = [M][u][Ar]"^ 

To impose boundary conditions in the r and 6 directions, the last and first rows of the 
operator matrices [A\ and [B] must be modified to include the known boundary values 
on the right-hand side of the equation ( Trefethen 2001| ) . Dirichlet boundary conditions 
are applied to the Helmholtz equation for v„ (see Section 3.5.21. Neumann boundary 
conditions are applied to the Poisson equation for pn^s (see Section C.2|. 
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